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Executive  Summary 


The  accurate  simulation  and  modeling  of  massively  separated  flows  comprises  one  of  the  primary 
obstacles  in  application  of  Computational  Fluid  Dynamics  (CFD)  as  a  more  routinely  used  tool  in 
engineering  design  and  analysis.  The  work  reported  here  represents  the  first  step  in  a  program  in 
which  the  long-term  objective  is  the  development  of  a  CFD  tool  for  predicting  aircraft  spin.  Such 
a  tool  will  offer  numerous  advantages  over  current  approaches  to  predicting  spin  characteristics, 
as  well  as  other  related  phenomena.  The  particular  focus  is  on  modeling  massive  separations  us¬ 
ing  Detached-Eddy  Simulation  (DES),  a  hybrid  method  that  combines  Reynolds-averaged  Navier- 
Stokes  (RANS)  approaches  with  Large  Eddy  Simulation  (LES).  The  primary  vehicle  for  assessing 
the  simulation  methodology  was  prediction  of  the  flow  around  a  model  of  an  aircraft  forebody, 
in  crossflow  with  the  freestream  at  angle  of  attack.  The  flow  at  angle  of  attack  is  sensitive  to  the 
Reynolds  number  with  experimental  measurements  showing  that  the  side  force  reverses  sign  with 
increases  in  Reynolds  number,  a  reversal  being  crucial  as  it  propels  the  spin  of  an  actual  forebody. 

Side-force  reversal  is  dependent  on  the  state  of  the  separating  boundary  layer,  i.e.,  laminar  or 
turbulent,  and  in  the  first  phase  of  the  work,  treatment  of  boundary  layer  separation  is  assessed 
using  calculations  of  the  two-dimensional  separated  flows  around  two  configurations:  a  rounded- 
comer  square  with  the  corner  radius  of  the  model  equal  to  1/4  of  the  body  width/diameter,  and 
a  2:1  ellipse.  The  incompressible  flow  around  the  configurations  is  predicted  using  a  fractional 
step  method  developed  for  application  in  curvilinear  coordinates.  A  grid  generation  scheme  was 
developed  using  the  control  technique  of  Hsu  and  Lee  (1991)  for  creation  of  structured  meshes  and 
used  for  gridding  both  geometries.  This  method  assures  orthogonality  of  the  grid  at  boundaries 
and  control  over  the  spacing  to  the  first  mesh  point  within  the  domain.  The  simulation  approach 
for  this  phase  of  the  work  was  based  on  solution  of  the  unsteady  RANS  (URANS)  equations  with 
the  Reynolds  stress  closed  using  the  Spalart-Allmaras  one-equation  model.  Characteristics  of  the 
separating  boundary  layer  are  established  by  the  inlet  conditions.  For  flows  with  turbulent  bound¬ 
ary  layer  separation,  a  small  level  of  eddy  viscosity  is  prescribed  at  the  inlet  to  the  computational 
domain,  sufficient  to  activate  the  turbulence  model  as  the  fluid  contacts  the  forebody.  Separation 
of  laminar  boundary  layers  is  accomplished  using  the  “tripless”  approach  of  Travin  et  al  (2001)  in 
which  the  initial  eddy  viscosity  within  the  domain  is  non-zero  and  the  inlet  eddy  viscosity  is  zero. 
The  separating  boundary  layer  is  laminar,  recirculation  in  the  wake  of  the  body  is  sufficient  to 
sustain  the  model.  Simulation  results  for  each  configuration  show  that  both  approaches  are  viable, 
with  the  distinctly  different  boundary  layer  separation  characteristics  achieved,  leading  to  changes 
in  the  pressure  coefficient  and  skin  friction  distributions  around  the  body. 

In  the  second  phase  of  the  work,  URANS  and  DES  are  used  to  predict  the  flow  around  the 
forebody  cross-section  modeled  by  the  rounded-comer  square,  the  same  cross-section  considered 
in  the  first  phase  of  the  study.  The  inlet  velocity  is  inclined  at  10°  to  the  main  flow,  the  configu- 
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ration  modeling  the  massively  separated  flow  around  the  forebody  of  a  jet  fighter  rotating  at  high 
angle  of  attack.  The  geometry  is  uniform  (extruded)  along  the  statistically  homogeneous  spanwise 
coordinate  for  which  periodic  boundary  conditions  were  applied.  Simulations  are  performed  at  a 
sub-critical  Reynolds  number  of  105  (based  on  the  freestream  speed  and  body  width/diameter)  for 
which  the  separating  boundary  layer  is  laminar  and  a  super-critical  Reynolds  number  of  8  x  105  for 
which  the  separating  boundary  layer  is  turbulent.  Between  these  Reynolds  numbers  experimental 
measurements  show  a  reversal  of  the  side  force.  DES  predictions  are  evaluated  using  experimen¬ 
tal  measurements  and  contrasted  with  the  URANS  predictions.  The  role  of  the  grid  is  assessed, 
initially  through  grid  refinement  performed  using  structured  grids.  Subsequently,  the  use  of  un¬ 
structured  grids  for  simulations  of  unsteady,  eddy -resolving  turbulent  flows  is  also  assessed.  The 
unstructured  grids  are  comprised  of  nearly  isotropic  tetrahedra  away  from  solid  surfaces  while  the 
boundary  layers  are  comprised  of  prisms.  DES  predictions  show  that  following  flow  detachment, 
a  chaotic  and  three-dimensional  wake  rapidly  develops.  The  temporal  evolution  of  the  streamwise 
and  lateral  (side)  forces  acting  on  the  body  exhibit  strong  modulation  due  to  the  spanwise  variation 
of  the  flow.  Grid  refinement  deepens  the  structure  of  the  resolved  range  of  turbulent  scales,  the 
predictions  on  unstructured  meshes  are  also  demonstrated  to  be  equally  accurate  as  those  on  struc¬ 
tured  grids.  For  the  super-critical  flow,  the  pressure  distribution  is  close  to  the  measured  values, 
both  the  streamwise  and  side  forces  are  in  adequate  agreement  with  measurements,  and  the  effect 
of  numerical  parameters  are  well-understood.  For  the  sub-critical  flow,  DES  side-force  predictions 
do  not  follow  the  experimental  measurements  far  enough  to  achieve  reversal.  Possible  causes  for 
the  discrepancy  are  discussed. 
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1  Introduction  and  Overview 


Knowledge  of  the  spin  and  recovery  characteristics  of  modem  aircraft  is  crucial  at  a  variety  of 
levels,  including  maneuverability,  control  strategies,  and  ultimately  design.  One  of  the  most  sig¬ 
nificant  factors  affecting  spin  characteristics  for  modem  fighters  is  the  forebody,  with  its  complex 
vortical  flows  and  long  moment  arm.  Laboratory  measurements  of  spin  characteristics  are  of  lim¬ 
ited  utility  since  it  is  not  possible  to  resolve  important  Reynolds  number  effects  because  of  the 
range  of  available  tunnels.  Numerical  simulation,  therefore,  provides  an  important  tool  that  should 
ultimately  provide  higher-fidelity  evaluations  of  aircraft  spin  than  current  approaches. 

Predicting  the  physics  of  spin  has  remained  a  substantial  challenge  to  modeling  approaches, 
e.g.,  the  flows  are  characterized  by  unsteady  massive  separation  and  the  dependence  on  transition. 
Simulation  techniques  used  for  prediction  of  these  phenomena  at  high  Reynolds  numbers  have 
traditionally  relied  on  Reynolds-averaged  approaches  which  are  unable  to  represent  these  physics 
to  sufficient  accuracy.  The  enormous  computational  requirements  to  resolve  boundary  layers  using 
techniques  which  incorporate  more  flow  details,  such  as  Large  Eddy  Simulation  (LES),  will  prevent 
their  application  to  whole  domains  at  the  Reynolds  number  ranges  encountered  in  applications  for 
the  foreseeable  future  (e.g.,  see  Spalart  et  al.  1997,  Spalart  2000). 

Unfortunately,  the  performance  of  numerical  models  has  been  inaccurate  in  most  instances 
due  to  their  inability  to  accurately  predict  the  complex  and  unsteady  effects  associated  with  spin. 
Vortical  flows,  crossflow  separations,  and  sensitivity  of  forces  and  moments  to  Reynolds  number 
greatly  challenge  modeling  approaches.  These  factors  also  supply  the  overall  motivation  for  the 
present  fundamental  investigations  and  the  need  to  develop  and  assess  improved  techniques  for 
predicting  complex,  separated  flows  at  high  Reynolds  numbers. 

1.1  Background:  Spin 

The  characteristic  motion  of  an  aircraft  following  stall  in  which  the  plane  descends  rapidly  along 
a  more  or  less  helical  trajectory  is  known  as  the  spin.  Spin  has  been  an  important  problem  for 
civil  and  military  aircraft  since  the  beginning  of  flight  (see  Durand  1934  and  Anderson  1979  for 
historical  reviews).  While  important,  prediction  of  spin  characteristics  is  difficult  and  is  dependent 
upon  a  large  number  of  aerodynamic  parameters,  including  slope  of  the  lift  curve  top,  rolling 
moment,  roll  and  yaw  damping,  pitching  and  yawing  moments,  etc.  While  stable  at  low  angle  of 
attack,  some  of  these  rotary  derivatives  can  change  sign  at  stall,  remain  unstable  at  high  attack 
angles  and  promote  spin. 

Early  spin  models  were  developed  from  force  and  moment  equilibrium  under  steady  conditions 
(e.g.,  see  Tischler  &  Barlow  1980  and  references  therein).  The  objective  was  to  identify  those 
factors  most  affecting  the  balance  of  forces  and  moments  during  a  steady  spin  and  to  use  linear 
analyses  to  consider  dynamic  behavior  for  states  near  the  steady  spin.  Linear  theories,  while  able 
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to  predict  some  aspects  of  steady  spins,  cannot  account  for  important,  and  highly  nonlinear,  effects 
such  as  sideslip. 

Most  of  what  is  known  about  spin  modes  has  been  derived  from  rotary  balance  data  (e.g.,  see 
Bihrle  et  al  1978,  Hultberg  et  al  1980).  Recent  measurements  and  analysis  of  rotary  balance 
data  obtained  at  the  Air  Force  Research  Laboratory  (AFRL),  while  not  directly  considering  spin, 
provides  an  illustration  of  the  complex  fluid  mechanics  which  challenge  predictions  (Jenkins  et  al 
1996,  Jobe  et  al  1996,  Grismer  &  Jenkins  1997,  Jenkins  1997).  Reported  in  these  investigations 
are  in-depth  studies  of  a  65  degree  delta  wing  at  high  angles  of  attack.  The  main  focus  was 
on  the  complex  development  of  the  leading-edge  vortices  in  regimes  where  vortex  breakdown 
occurs.  Aerodynamic  reactions  are  sensitive  to  breakdown  and  related  phenomena  which  occur 
over  disparate  time  scales.  Because  linear  models  cannot  account  for  these  effects,  subsequent 
efforts  have  focused  on  development  of  nonlinear  models  with  time  constants  in  aerodynamic 
response  determined  from  measurements  (e.g.,  see  Tobak  &  Chapman  1985,  Troung  &  Tobak 
1990,  Myatt  1997). 

Simulation  work  relevant  to  the  the  research  reported  here  includes  the  solutions  of  the  un¬ 
steady,  three-dimensional  compressible  Navier-Stokes  equations  of  the  flow  around  the  delta  wing, 
but  at  lower  Reynolds  number,  than  considered  in  the  AFRL  rotary  measurements.  The  objective 
was  to  determine  if  yaw  rate  effects  evident  in  coning  motions  from  rotary  balance  tests  will  sig¬ 
nificantly  change  the  mean  location  of  vortex  breakdown.  Tromp  (1998)  investigated  effects  of 
roll  and  yaw  rates  by  considering  the  responses  from  a  coning  motion  and  a  pure  roll-rate  helical 
motion.  Tromp  found  that  the  movement  in  vortex  breakdown  location  with  yaw  rate  was  rela¬ 
tively  small,  suggesting  that  using  coning  motion  to  estimate  the  rate  effects  on  critical  states  may 
provide  reasonable  estimations. 

An  effect  recognized  to  be  important,  but  not  considered  in  the  work  summarized  above  is 
that  of  the  Reynolds  number.  Though  spin  characteristics  are  known  to  exhibit  strong  sensitivity 
to  Reynolds  number  (e.g.,  see  McCormick  1981),  there  have  been  relatively  few  investigations 
of  Reynolds  number  effects  in  spin  or  spin-related  flows.  One  attempt  is  the  work  reported  by 
Fritz  (1995),  who  used  RANS  computations  to  predict  spin  characteristics  of  generic  forebody 
cross-sections.  Two-dimensional  solutions  were  compared  to  two-dimensional  wind  tunnel  data  of 
Polhamus  et  al  (1959)  and  a  strip  theory  method  was  used  to  integrate  the  solutions  to  compare 
with  data  gathered  for  four  generic,  rotary  balance  models.  Fritz  (1995),  while  able  to  predict  the 
proper  spin  in  super-critical  Reynolds  number  regimes,  was  unable  to  predict  the  proper  spin  in 
sub-critical  regimes.  The  major  difficulties  encountered  were  correctly  and  consistently  predicting 
the  location  of  flow  separation  and  the  inability  of  turbulence  models  to  handle  the  varying  length 
scales  of  separated,  bluff  body  flows. 

The  lack  of  reliable  methods  for  handling  Reynolds  number  variations  severely  limits  CFD 
for  predicting  spin  and  improved  approaches  have  long  been  sought.  Experimental  resolution  of 
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Reynolds  number  effects  requires  costly  and  time-consuming  tests,  probably  using  more  than  one 
wind  tunnel.  Models  of  aerodynamic  response  require  some  calibration  using  either  experiments  or 
simulations.  Unfortunately,  previous  computational  efforts  have  not  considered  Reynolds  number 
effects  or  have  been  unable  to  predict  spin  characteristics  over  a  range  of  Reynolds  numbers  using 
available  strategies,  i.e.,  RANS  techniques  without  clear  control  of  transition.  The  drawback  of 
current  simulation  strategies  based  on  solution  of  the  Reynolds-averaged  Navier-Stokes  equations 
is  that  accurate  prediction  of  spin  requires  accurate  prediction  of  geometry-dependent,  unsteady 
three-dimensional  turbulent  motions  which  are  very  important  in  the  separated  flows  characteriz¬ 
ing  spin.  These  eddies,  arguably,  are  what  defeats  traditional  RANS  turbulence  models,  of  any 
complexity.  Presented  next  is  an  overview  of  methods  used  for  predicting  complex  flows  at  high 
Reynolds  numbers,  an  overview  that  is  used  to  provide  the  context  for  the  technique  forming  the 
backbone  for  this  research:  Detached-Eddy  Simulation. 

1.2  Background:  Hybrid  methods  for  predicting  high  Reynolds  number 
flows 

Most  high-Reynolds  number  predictions  are  currently  obtained  from  solutions  of  the  Reynolds- 
averaged  Navier-Stokes  (RANS)  equations.  While  the  most  popular  RANS  models  appear  to  yield 
predictions  of  useful  accuracy  in  attached  flows  as  well  as  some  with  shallow  separations,  RANS 
predictions  of  massive  separations  have  typically  been  unreliable.  RANS  models,  calibrated  in 
thin  shear  layers,  appear  unable  to  consistently  represent  to  sufficient  accuracy  the  geometry- 
dependent,  chaotic  and  unsteady  features  of  massively  separated  flows.  This  is  the  case  even 
with  unsteady  RANS,  which  is  often  capable  of  producing  vortex  shedding,  though  shedding  from 
URANS  predictions  is  typically  non-chaotic  and  exaggerated  in  amplitude  (e.g.,  see  Travin  et  al 
2001).  URANS  predictions  of  drag,  for  example,  have  typically  been  too  high,  though  it  should  be 
recognized  that  this  pessimistic  view  of  URANS  is  not  a  universal  one  (Durbin  1995). 

The  relatively  poor  performance  of  RANS  models  has  motivated  the  increased  application  of 
Large  Eddy  Simulation  (LES).  Away  from  solid  surfaces,  LES  is  a  powerful  approach,  providing 
a  description  of  the  large,  energy-containing  scales  of  motion  that  are  typically  dependent  on  ge¬ 
ometry  and  boundary  conditions.  When  applied  to  boundary  layers,  however,  the  computational 
cost  of  whole-domain  LES  does  not  differ  significantly  from  that  of  Direct  Numerical  Simulation 
(DNS)  (Spalart  et  al  2000).  The  “large  eddies”  close  to  the  wall  are  physically  small  in  scale,  in  a 
whole-domain  LES  (i.e.,  without  wall  modeling),  these  small  near-wall  structures  are  dynamically 
important  and  must  be  resolved  in  order  to  accurate  predict  boundary  layer  evolution.  In  boundary 
layers  without  sufficient  resolution  of  the  wall-layer  structures,  LES  predictions  of  boundary  layer 
growth  and/or  separation  will  be  adversely  affected. 

These  and  other  considerations  provided  the  motivation  for  development  of  Detached-Eddy 
Simulation  (DES)  by  Spalart  et  al  (1997).  DES  belongs  to  the  class  of  hybrid  methods  which 
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attempt  to  combine  RANS  and  LES,  exploiting  the  accuracy  and  efficiency  of  RANS  approaches  in 
the  thin  shear  layers  not  far  from  the  calibration  range  of  the  model,  and  LES  for  direct  resolution  of 
the  geometry-dependent  and  three-dimensional  eddies  in  regions  where  such  a  treatment  is  desired. 
DES  was  developed  independently  of  a  different  approach  proposed  by  Speziale  (1998)  also  aimed 
at  combining  RANS  and  LES.  Speziale  (1998)  proposed  damping  the  Reynolds  stress  predicted 
by  a  second-moment  closure  approach,  (uiUj)^M\  by  a  factor  a, 

(■ UiUj )  =  a(v,iUj)(M ^ ,  (1) 

where  a  <  1.  The  Reynolds  stress  ( u^Uj )  on  the  left-hand  side  (1)  is  then  supplied  to  the  momen¬ 
tum  equations.  Speziale  (1998)  proposed  that  the  damping  a  be  determined  based  on  a  measure  of 
the  local  resolution  compared  to  the  Kolmogorov  scale. 


a  = 


1  -  qxv(-P(L£,/Lk) 


n 


(2) 


where  LA  is  a  lengthscale  related  to  the  grid  size,  Lk  is  the  Kolmogorov  lengthscale  and  /?  and 
n  are  adjustable  parameters.  Two  limits  of  (2)  are  well  defined.  The  first  limit  corresponds  to 
the  “RANS  limit”  of  the  model  for  which  LA  >>  Lk  in  which  case  a  — >  1  and  the  Reynolds 
stress  supplied  to  the  momentum  equations  is  identical  to  that  obtained  from  the  particular  second- 
moment  closure  being  employed.  The  second  limit  corresponds  to  the  “DNS  limit”  for  which 
LA  <<  Lk  in  which  a  ->  0  and  the  Reynolds  stress  supplied  to  the  governing  equations  vanishes. 
Intermediate  between  these  regimes  the  damping  factor  a  will  be  less  than  unity,  corresponding  to 
only  a  fraction  of  the  model  stress  a(uiUj)^  predicted  using  the  second-moment  closure  being 
supplied  to  the  Navier-Stokes  equations.  In  this  intermediate  range,  therefore,  part  of  the  turbulent 
stress  is  modeled,  the  remaining  part  to  be  resolved  directly  on  the  mesh.  This  intermediate  regime 
corresponding  to  the  “LES  regime”  of  the  model.  An  important  question  is  the  response  of  the 
model  in  its  “LES  regime”  for  computations  in  which  LA  >>  Lk,  the  criteria  also  used  to  define 
the  “RANS  region”,  i.e.,  the  form  of  the  model  should  be  fundamentally  different  in  the  RANS  and 
LES  limits  though  it  appears  from  (1)  and  (2)  that  there  could  be  some  ambiguity  when  LA  >> 
Lk-  Von  Terzi  and  Fasel  (2002)  have  adopted  the  basic  outline  of  the  above  approach,  one  of  the 
main  modifications  being  the  introduction  of  a  “contribution  function”  which  plays  a  similar  role 
as  a  in  (2),  dictating  the  level  of  the  turbulent  stress  to  be  modeled  versus  resolved. 

Batten  et  al  (2000)  have  developed  a  hybrid  RANS-LES  method  referred  to  as  Limited  Nu¬ 
merical  Scales  (LNS).  The  approach  appears  to  be  very  similar  to  DES,  the  method  reducing  to  a 
RANS  treatment  of  the  boundary  layer  and  LES  away  from  solid  surfaces.  The  relation  (1)  com¬ 
prised  the  basis  for  the  initial  form  of  the  model,  though  rather  than  defining  a  using  (2),  Batten  et 
al  (2000)  propose  to  base  a  essentially  on  a  ratio  of  eddy  viscosities  in  the  RANS  and  LES  regions, 

_  min ({Lv)LESi  {Lv)ea) 
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where  ( Lv)les  represents  a  lengthscale-velocity  product  obtained  from  an  LES  model  and  (Lv)ba 
represents  a  lengthscale-velocity  product  obtained  from  an  ensemble-averaged  model. 

As  described  in  greater  detail  in  §3.1,  DES  has  RANS  behavior  near  the  wall  and  becomes 
a  Large  Eddy  Simulation  in  the  regions  away  from  solid  surfaces  provided  the  grid  density  is 
sufficient.  The  formulation  is  based  on  a  modification  to  the  Spalart-Allmaras  one-equation  model 
(Spalart  and  Allmaras  1994,  referred  to  as  S-A  throughout).  Some  of  the  principle  advantages  of 
DES  are  that  the  technique  is  non-zonal  and  computationally  feasible  for  high  Reynolds  number 
prediction,  but  also  resolves  time-dependent,  three-dimensional  turbulent  motions  as  in  LES. 

A  recent  review  of  some  of  the  applications  of  DES  can  found  in  Strelets  (2001).  DES  has  pre¬ 
dicted  to  much  higher  accuracy  than  traditional  approaches  the  flow  around  an  airfoil  at  high  angle 
of  attack,  another  necessary  pre-requisite  for  prediction  of  spin  (Shur  et  cil.  1999).  DES  has  also 
been  used  to  successfully  predict  the  flow  around  a  sphere,  an  important  benchmark,  resolving  the 
time- dependent  forces  around  the  body  to  a  higher  accuracy  than  methods  based  on  the  Reynolds- 
averaged  Navier-Stokes  (RANS)  equations,  and  in  a  Reynolds  number  range  higher  than  can  be 
realistically  attempted  using  full-domain  LES  (Constantinescu  et  cil.  2002).  An  important  outcome 
of  these  studies  is  that  the  cost  of  DES  exhibits  a  weak  dependence  on  Reynolds  number,  no  worse 
than  RANS  methods. 

Because  the  boundary  layer  is  predicted  by  the  S-A  model,  the  empirical  input  to  the  technique, 
as  true  for  all  hybrid  methods,  is  not  small  and  in  DES,  for  example,  the  prediction  of  boundary 
layer  separation  is  under  control  of  the  RANS  model.  Assessing  the  technique,  therefore,  remains 
an  important  task,  not  only  with  respect  to  aspects  such  as  predicting  boundary  layer  separation 
(the  primary  consideration  of  the  first  part  of  this  work),  but  also  the  response  of  the  method  to  the 
numerical  approach,  e.g.,  the  grid  and  timestep,  which  are  some  of  the  other  aspects  scrutinized  in 
the  second  part  of  the  work. 

2  Objectives  and  summary  of  the  approach 

The  research  described  in  this  report  is  a  pre-cursor  to  the  ultimate  application  which  is  prediction 
of  the  spin  characteristics  of  full  aircraft  at  flight  Reynolds  numbers.  The  flow  fields  characterizing 
spinning  aircraft  are  massively  separated,  providing  a  “natural”  application  for  DES  and  therefore 
the  principle  objective  was  to  apply  and  assess  the  method  in  a  flow  important  to  understanding  and 
accurately  predicting  spin.  Though  a  natural  application  for  the  model,  calculations  of  complex 
configurations  at  high  Reynolds  numbers  challenge  the  entire  computational  approach,  including 
numerical  aspects  related  to  features  such  as  grid  generation  and  grid  type  (i.e.,  structured  or  un¬ 
structured),  to  aspects  impacting  the  physical  modeling,  an  issue  of  primary  importance  being  the 
treatment  of  the  separating  boundary  layer. 

The  work  is  divided  into  two  parts.  For  spin  prediction,  a  key  issue  is  the  sign  of  the  side  force, 
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and  whether  the  side  force  is  spin-damping  or  spin-propelling.  The  side  force,  both  in  magni¬ 
tude  and  sign,  is  dependent  the  state  of  the  separating  boundary  layer,  i.e.,  whether  the  separating 
boundary  layer  is  laminar  or  turbulent.  In  the  context  of  hybrid  methods,  for  which  boundary  layer 
separation  is  under  control  of  the  RANS  model,  the  type  of  transition  is  controlled  by  the  eddy  vis¬ 
cosity  using  the  S-A  model.  In  the  first  phase  of  the  work,  calculations  of  two-dimensional  flows 
were  used  to  investigate  approaches  for  controlling  the  type  of  boundary  layer  separation  and  to 
assess  viability  for  the  more  computationally-intensive  three-dimensional  calculations  performed 
in  the  second  phase  of  the  work. 

In  the  second  phase  of  the  work  the  factors  included  in  the  investigation  include  the  role  of 
the  grid,  assessment  of  the  flow  solver  and  treatment  of  boundary  layer  detachment  for  three- 
dimensional  configurations.  Given  the  end  application  of  full-aircraft  configurations,  unstructured 
grids  form  an  integral  component  of  the  present  approach  and  therefore  an  unstructured-grid  solver 
was  employed  for  numerical  solution  of  the  Navier-Stokes  equations  for  the  three-dimensional 
calculations. 

The  disadvantage  of  the  current  unstructured  approach  is  that  the  numerical  procedure  is  only 
second-order  accurate  in  time  and  space  and  stabilized  via  non-linear  (TVD)  numerical  dissipa¬ 
tion.  Related  investigations  have  shown  that  the  artificial  dissipation  associated  with  the  numerical 
scheme  can  be  as  large  as  that  represented  by  the  turbulence  model  and  therefore  care  must  be 
exercised  in  application  of  these  methods  to  eddy-resolving  simulations  such  as  LES  (e.g.,  see 
Mittal  and  Moin  1997).  This  in  turn  motivates  another  goal  of  the  present  effort  -  to  explore  the 
accuracy  of  the  current  second-order  method  on  both  structured  and  unstructured  grids  for  DES 
applications.  This  is  not  a  simple  matter  of  verifying  the  order  of  accuracy,  which  is  difficult  to 
define  and  predict  in  LES  and  especially  hybrid  methods.  Nevertheless,  the  primary  tool  for  such 
a  study  remains  grid  refinement. 

3  Approach 

3.1  Detached-Eddy  Simulation 

Detached  Eddy  Simulation  (DES)  is  a  hybrid  technique  first  proposed  by  Spalart  et  cil  (1997)  for 
prediction  of  turbulent  flows  at  high  Reynolds  numbers  (see  also  Spalart  2000).  The  motivation  for 
developing  DES  was  to  combine  the  most  favorable  aspects  of  RANS  and  LES,  i.e.,  the  acceptable 
predictions  using  RANS  models  of  thin  shear  layers  and  LES  for  resolution  of  time-dependent, 
three-dimensional  large  eddies  which  are  typically  geometry-dependent.  RANS  models  are  less 
accurate  in  separated  regions  than  in  the  thin  shear  layers  where  they  are  calibrated.  In  these 
regions  LES  is  very  attractive  since  the  large  scales  can  be  resolved  without  the  vast  increases  in 
grid  resolution  necessary  in  LES  of  boundary  layers. 

The  DES  formulation  is  based  on  a  modification  to  the  S-A  model  such  that  it  reduces  to 
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RANS  close  to  solid  surfaces  and  to  LES  away  from  the  wall  (Spalart  et  al.  1997).  The  S-A  RANS 
model  of  Spalart  and  Allmaras  (1994)  is  summarized  below  along  with  a  discussion  of  the  DES 
formulation.  Additional  discussion  can  be  found  in  Spalart  (2000)  and  Strelets  (2001). 

In  the  S-A  RANS  model,  a  transport  equation  is  used  to  compute  a  working  variable  used  to 
form  the  turbulent  eddy  viscosity, 


Du 

Dt 


—  CblS  V  CVJ  i  f w  f t2 J  ^ 


1 


+  —  [V  •  {(v  +  u)Vu)  +  C(,2  (Vi7)  j  +  fa  AU 


where  v  is  the  working  variable.  The  eddy  viscosity  vt  is  obtained  from, 

'3  _  v 

X  v  > 


Ui  —  v  fv  1,  ,fvl 


X3  +  c3i  ’ 


where  v  is  the  molecular  viscosity.  The  production  term  is  expressed  as, 


S  =  fv3S  +  -^pfv2  , 


(4) 

(5) 

(6) 


fv  2  = 


fv  3  ~ 


(1  +  x/di)(1  ~  fv 2) 

X 


(7) 


where  S  is  the  magnitude  of  the  vorticity.  The  production  term  as  written  in  (6)  differs  from 
that  developed  in  Spalart  and  Allmaras  (1994)  via  the  introduction  of  fv 3  and  re-definition  of  /„ 2- 
These  changes  do  not  alter  predictions  of  fully  turbulent  flows  and  have  the  advantage  that  for 
simulation  of  flows  with  laminar  separation,'  spurious  upstream  propagation  of  the  eddy  viscosity 


into  attached,  laminar  regions  is  prevented.  The  function  fw  is  given  by, 


fw  =  9 


1  +  c6 


w3 


9°  + c! 


w3 


1/6 


9  =  r  +  cw2  (r6  -r), 


r  = 


(8) 


The  function  fa  is  defined  as, 

ft 2  =  ci3exp(-ct4X2)  •  (9) 

The  trip  function  ftl  is  specified  in  terms  of  the  distance  dt  from  the  field  point  to  the  trip,  the  wall 
vorticity  cj;  at  the  trip,  and  AU  which  is  the  difference  between  the  velocity  at  the  field  point  and 
that  at  the  trip, 

fti  =  Cu.9iexP  (^~Cl2~£jp  [d2  +  g2d2]  ^ 

where  gt  =  min(0.1.  AU/utAx)  and  Ax  is  the  grid  spacing  along  the  wall  at  the  trip.  The  wall 
boundary  condition  is  V  =  0.  As  described  in  §4,  the  inlet  condition  on  v  was  prescribed  for  either 
fully  turbulent  boundary  layers  prior  to  separation  or  with  zero  upstream  values  of  v  that  yielded 
laminar  boundary  layers  prior  to  separation.  The  constants  are  c&i  =  0.1355,  a  —  2/3,  c& 2  =  0.622, 
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K  —  0.41,  Cwi  —  C-bi/ K?  +  (l  +  Cbz) / & i  C-vj2  —  0.3,  Cw 3  —  2,  Cv i  —  7.1,  Cv 2  5,  Ct\  1,  Ct2  2, 

c*3  =  1.1,  andct4  =  2. 

The  DES  formulation  is  obtained  by  replacing  the  distance  to  the  nearest  wall,  d,  by  d,  where 
d  is  defined  as, 

d  =  min(d,  Cdes&)  »  (H) 

with 

A  =  max(Ax,  A y,  Az) .  (12) 

where  Ax,  Ay,  and  Az  are  the  grid  spacings.  In  “natural”  applications  of  DES,  the  wall-parallel 
grid  spacings  (e.g.,  streamwise  and  spanwise)  are  at  least  on  the  order  of  the  boundary  layer  thick¬ 
ness  and  the  S-A  RANS  model  is  retained  throughout  the  boundary  layer,  i.e.,  d-d.  Conse¬ 
quently,  prediction  of  boundary  layer  separation  is  determined  in  the  ‘RANS  mode’  of  DES.  Away 
from  solid  boundaries,  the  closure  is  a  one-equation  model  for  the  SGS  eddy  viscosity.  When  the 
production  and  destruction  terms  of  the  model  are  balanced,  the  length  scale  d  —  Cqes&  in  the 
LES  region  yields  a  Smagorinsky  eddy  viscosity  ?oc5A2.  Analogous  to  classical  LES,  the  role  of 
A  is  to  allow  the  energy  cascade  down  to  the  grid  size;  roughly,  it  makes  the  pseudo-Kolmogorov 
length  scale,  based  on  the  eddy  viscosity,  proportional  to  the  grid  spacing.  “Wall  distance”  is  the 
distance  to  the  nearest  wall,  which  can  be  clearly  defined  even  in  complex,  fully  three-dimensional 
geometries  with  multiple  length  scales.  The  additional  model  constant  C des  =  0*65  was  set  in 
homogeneous  turbulence  (Shur  et  al.  1999)  and  is  used  without  modification  in  this  work. 

There  are  several  advantages  to  the  DES  formulation  described  above.  The  transition  between 
RANS  and  LES  is  seamless  in  a  formulation  sense:  single  equation,  no  explicit  declaration  of 
RANS  versus  LES  zones.  The  solutions  shown  below  show  that  the  DES  is  seamless  in  an  appli¬ 
cation  sense,  i.e.,  no  artificial  transitions  between  the  solution  domains.  There  is  one  solution  field, 
all  coupled  by  the  Navier-Stokes  equations,  and  the  solution  has  a  different  character  in  different 
regions. 

A  further  advantage  of  the  length  scale  switch  in  DES  is  that  it  allows  the  user  to  ‘steer’  the 
physics  where  needed,  i.e.,  finer  resolution  in  regions  where  the  additional  flow  detail  from  an  LES 
treatment  is  desired.  The  length  scale  switch  in  the  model  has  the  effect  of  drawing  down  the  eddy 
viscosity  so  the  modeled  equation  can  allow  development  of  instabilities,  which  then  feed  a  cascade 
to  smaller  scales,  limited  by  the  grid  as  in  canonical  LES.  Regions  of  the  flow  not  far  from  the  thin 
shear  layer  approximation  are  very  expensive  to  resolve  using  LES,  and  can  be  adequately  handled 
by  the  S-A  model  within  a  RANS  treatment.  Important  for  applications  is  that  the  increase  in  cost 
of  DES  with  Reynolds  number  is  weak,  similar  to  RANS,  unlike  in  full-domain  LES  in  which 
there  is  an  enormous  increase  in  cost  because  of  the  rapid  increase  in  resolution  requirements  near 
walls  (see  also  Nikitin  et  al  2000). 
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3.2  Numerical  methods 

3.2.1  Incompressible  flow  solver 


The  approach  to  investigation  of  boundary  layer  separation  characteristics  for  the  first  phase  of  the 
study  was  based  on  solution  of  the  incompressible  Reynolds-averaged  Navier-Stokes  equations.  A 
code  for  implementation  of  a  fractional-step  numerical  method  was  developed,  the  basic  aspects 
of  the  numerical  approach  are  outlined  in  this  section. 

The  unsteady  incompressible  Navier-Stokes  equation  written  in  tensor  form  are  expressed  as. 


duj 

dxi 


=  0, 


(13) 


dui  duiiij 

Ht  + 


1  dp  d 


+ 


dxj 


dui  ^  duj 


I dxj  dxi  j 


(14) 


*>%  U'j 

dxj  p  dxi 

where  the  velocity  is  denoted  u* ,  the  pressure  p,  density  p,  and  v  is  the  fluid  kinematic  viscosity. 
Non-dimensionalizing  the  above  using  a  characteristic  length  L  and  characteristic  velocity  Woo 
results  in  a  re-definition  of  the  lengths,  time,  and  dependent  variables, 


Ui  —  U-i  J  Uq 


Xi  =  Xi/L , 


P 


P 

Pulo 


t  -  tVtQQ  j  L  , 


and  the  corresponding  form  of  the  non-dimensionalized  equations, 


dui  d 
dt  +  dXi{UiUj) 


dp  d  v  dui  duj 

dxi  dxj  UoqL  i  dxj  dxi  J 


(15) 


(16) 


The  numerical  approach  to  approximation  of  the  above  system  is  based  on  integration  of  (16) 
over  an  arbitrary  (fixed)  control  volume  V , 


L>  -  Lb~^)iv + L-k[~v)iV 


+ 


!  — 

1 

dui 

duj 

jv  dxj 

Re 

dxj 

dxi 

(17) 


with  the  Reynolds  number  Re  =  u^L/v.  The  equations  can  be  cast  equivalently  in  vector  form 
as, 


[  dV  =  [  V  •  (-UU )dV  +  [  V  •  (-PI)dV 

Jv  dt  Jv  Jv 


+ 


[  V  ■  — 

Jv  ^ 


Vu  +  (Vuf 


dV, 


(18) 


where  uu  is  a  dyad  (second  order  tensor).  Using  Gauss’  theorem,  fv  V  •  AdV  =  Js  A  •  dS,  the 
equations  governing  the  constant-density  and  isothermal  flow  over  the  fixed  control  volume  V  with 
surface  S  can  finally  be  written  as, 

[  u  ■  dS  =  0  (19) 
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[  =  [  T  ■  dS 

Jv  9t  Js 


(20) 


where  T  is  given  in  operator  form  as, 


T  =  -uu  -  PI  +  ~  [Vu  +  (Vu)7]  , 


or  in  tensor  notation  expressed  as 


T  ij  —  UiUj  Pfiij  + 


dui  duj 
dxj  dxi 


(21) 


(22) 


Coordinate  transformation 

The  equations  (19)  and  (20)  are  solved  in  a  general  curvilinear  system  that  is  simplified  for  the 
present  investigations  to  two-dimensional  configurations,  i.e.,  with  a  uniform  (non-transformed) 
spanwise  coordinate  z.  The  first  step  in  solution  of  (19)  and  (20)  is  consideration  of  the  coordinate 
transformation  of  a  vector  R  from  the  Cartesian  system  to  a  curvilinear  non-orthogonal  system  as 

shown  schematically  in  Figure  la,  R(x,  y,  z)  ->  R(£,  V,  z). 

In  what  follows,  the  natural  local  base  of  the  transformed  curvilinear  coordinate  system  (£,  y,  z) 
is  denoted  eq  =  dr/dq,  q  =  r?,  2.  The  transformation  between  the  two  systems  is  achieved  by 

defining  a  new  base  in  terms  of  the  vector  product  of  the  natural  local  basis  eq ,  Sq  =  e9+1  x  eq+2 
(q,  cyclic  permutation).  Note  that  the  vector  S9  has  a  magnitude  equal  to  the  area  of  the  paral¬ 
lelogram  spanned  by  natural  local  basis  (e9+i,  e9+ 2)  and  is  perpendicular  to  the  plane  defined  by 
(e9+i,  e,i+2)-  Using  the  orthogonality  between  the  new  basis  Sq  and  its  reciprocal  base  Sg,  the  se¬ 
ries  operations  in  generalized  curvilinear  systems  as  outlined  below  will  naturally  yield  dependent 
variables  that  are  the  volume  fluxes,  rather  than  the  primitive  variables. 

At  an  arbitrary  point  (z,  j,  k)  within  the  domain,  the  area  tensor  is  defined  as, 

S  =  (S«,S^,SZ),  (23) 


where  the  vector  components  are  defined  as, 


dr  dr 

s 


__  dr  dr 

s  ~d~zxdV 


dr  dr 

S  ~dl^dy' 


or  in  a  more  compact  form, 


Sq  = 


dr 


dr 


(24) 


(25) 


d(q  +  l)  d(q  +  2) 

where  q  =  f ,  77, 2  are  in  cyclic  order.  The  vector  quantity  S9  has  the  magnitude  of  the  area  of  the 
face  and  direction  given  by  the  normal  of  the  face  (c.f..  Figure  lb). 

The  partial  derivative  of  the  position  vector  r(x(£,  ri),y(£,  1?) ,  2)  with  respect  to  £,  r/,  z  are  eval¬ 
uated  using  a  central  finite  difference  approximation.  For  example,  with  reference  to  Figure  lb, 


dr 


r[z(£,  77, 2),  y(£,  17),  z]e  -  [x(£,  77,  z),y{^y),  2] 
A£ 


(26) 
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The  area  vectors  S’  are  required  to  satisfy  geometric  constraints,  the  first  of  which  is  that  for 
any  computational  cell  within  the  domain  the  cell  should  be  closed,  i.e., 

J /dS  =  0,  (27) 


or  in  discrete  form, 


Vs’  =  o. 


(28) 


Q 

To  satisfy  (26)  to  machine  accuracy  in  the  computation,  dr/ dq  needs  to  be  consistently  evaluated 
with  respect  to  the  computational  grid  points.  In  the  present  computations,  the  area  vectors  Sq  and 
their  reciprocal  bases  evaluated  and  stored  at  cell  centers.  Area  vectors  at  cell  faces  are  obtained 
via  simple  arithmetic  averaging. 

The  volume  of  the  computational  cell  is  evaluated  from 


r  sj,  +  s?  +  s ; 


(r/  ne  rbsw)  j 


(29) 


with  the  geometric  constraint  on  the  volumes  specified  as, 


Y,V  =  V^nain  •  (30) 

cells 

Discretization  of  the  continuity  and  momentum  equations 

The  continuity  constraint  is  discretized  in  its  integral  form,  applied  to  a  given  computational 
volume  and  with  reference  to  Figure  lb  takes  the  form, 

(S«  •  u)e  -  (S«  •  u)w  +  (S’  •  u)n  -  (S’  •  u),  +  (S*  •  U )/  -  (S*  •  u)b  =  0 .  (31) 


The  form  (31)  suggests  dependent  variables  different  than  the  primitive  variables.  The  quantities 
in  (31)  represent  volume  fluxes  and  can  be  formally  defined, 

■ut  =  Sf  •  u  =  S^Ui  +  S^Uj  +  S^k  (32) 

uv  =  S’  •  u  =  5’u.j  +  SqUj  +  5’ti*;  (33) 

uz  =  Sz  •  u  =  SzUi  +  SzUj  +  Sz2uk  (34) 

where  v£,un,uz  are  the  volume  fluxes  over  the  rj,  z  faces  of  a  primary  control  volume,  and 

uu  Uj,  Uk  are  the  Cartesian  velocity  components. 

The  introduction  of  the  volume  fluxes  leads  to  the  continuity  constraint  expressed  as, 


4  +  K  -  v?s  +  u)  -  uzb  =  0 


(35) 
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where  uq  =  (u^,un.,uz).  An  advantage  of  the  use  of  volume  fluxes  and  continuity  constraint 
written  as  (35)  is  that  global  mass  conservation  is  easily  expressed  via, 


£  “'  =  0 

total  cell  faces 


(36) 


The  integral  form  of  the  momentum  equation  with  Cartesian  velocity  components  as  unknown 
variables  is  rewritten  below, 


IILPv  =  ILTdS- 

Approximating  the  above  over  an  arbitrary  (and  fixed)  control  volume  V  yields, 

V^-  =  VS?-T  =  F 

dt  ^ 


(37) 


(38) 


where  F  represents  the  total  flux  through  the  computational  cell. 

Using  the  relations  (32)-(34),  the  momentum  equations  along  direction  q  with  u^,u'1,uz  as 

variables  are,  _ 

S9  ■  (U^)  =  V ^  =  S9  •  F  (39) 

dt  at 

where  q  ==  £  or  77  or  z.  The  momentum  equations  in  the  form  (39)  are  those  considered  next  for 
numerical  approximation. 

Fractional  step  method 

Writing  (39  for  the  £  direction,  with  reference  to  Figure  lb, 


vM  =  S*-F  =  St'(£  S9-T]  =  LC, 


dt 

The  linear  operator  L %  is  split  into  several  parts, 

+  D£iex  , 


(40) 


(41) 


where  H %  are  the  nonlinear  convection  terms,  the  pressure  term,  D %  are  the  diffusion  terms  that 
are  advanced  implicitly,  and  D^ex  represents  the  diffusion  terms  advanced  explicitly.  The  decom¬ 
position  of  the  entire  diffusion  term  into  two  components  is  to  more  efficiently  facilitate  imple¬ 
mentation  of  the  fractional  step  method.  Essentially,  contains  terms  similar  to  d2u^/d£>2  + 

d2ut /drf  +  d2ut / dz 2,  while  D^ex  contains  the  rest  of  the  diffusion  terms,  those  depending  on  the 
cross  derivatives. 

The  explicit  portion  of  the  time  advance  is  achieved  using  second-order  Adams-Bashforth  for 
the  convection  terms, 

=  lH((uq)n  -  ,  (42) 
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the  pressure  term  advanced  using  implicit  Euler, 


R^P)  =  R({P)n+1 ,  (43) 

and  Crank-Nicholson  for  the  implicit  diffusion  terms, 

D^v?)  =  ^D?(^)n+1  +  .  (44) 

and  second-order  explicit  Adams -Bashforth  scheme  for  the  explicit  diffusion  terms, 

D^q)  =  lDi>ex(uT  -  Id^uT'1  (45) 

Finally,  the  temporal  discretization  of  (40)  can  be  written  for  the  three  components  of  the  momen¬ 
tum  equation  as, 

(^)"+1  ~  (u^P  =  _  ^(U9)n-lj  +  ^(pn+1)  (46) 

l\l  Z 

+i[3D?>ei(u«)"  -  D^uT-1}  +  \mn(r+1 + , 

1[3 Hv(u9)n  -  Hr,(u9)n~1}  +  Rn(Pn+l)  (47) 

+  i[3A),es(uT  -  A.esKr1]  +  \muT+l  +  Dv(uT]  , 

K(uZ)n+1A~------  =  l[3Hz(ull)n  -  Rz(Pn+1)  (48) 

ZA  6  Z 

+i[3 M«T  -  D,,a(uT-1)  +  i[A(»‘)"+1  +  D*(uT]  ■ 

Summary  of  the  fractional  step  method 

The  system  (31)  and  (46)-(48)  is  solved  using  a  fractional  step  method,  comprised  of  three 
sequential  operations:  computation  of  a  provisional  velocity  field  using  the  non-linear  and  viscous 
terms;  calculation  of  the  pressure  field  by  solving  the  Poisson  equation,  and  finally  projection  of 
the  intermediate  velocity  field  onto  a  divergence-free  space  at  the  new  time  step  using  the  pressure 
gradient.  Written  in  a  compact  form,  the  series  of  operations  are 


Au 

=  rhs  +  mbc 

(49) 

At  DG4>n+l 

=  Du  -  cbc 

(50) 

un+1 

=  u-A  tG<t>n+1 

(51) 

where  the  G  submatrix  is  the  discrete  gradient  operator  and  D  is  the  discrete  divergence  oper¬ 
ator.  The  exact  expressions  for  the  coefficient  submatrix  A  and  right-hand-side  submatrix  r  are 
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dependent  on  the  specific  temporal  and  spatial  discretizations.  The  unknown  discrete  velocity  and 
pressure  vectors  are  denoted  un+1  andpn+1.  Boundary  condition  vectors  mbc  and  cbc  are  defined 
such  that  when  combined  with  operators  G  and  D  the  original  system  with  appropriate  boundary 
conditions  can  be  fully  recovered. 

Given  the  transformation  to  a  general  curvilinear  system  as  outlined  above,  the  fractional  step 
method  as  implemented  takes  the  form  for  the  volume  fluxes  as, 

Afiq  =  (S?,  S’7,  S2)t  •  rhs  +  (Sf,  S’7,  SZ)T  •  mbc  (52) 

(Si,Sv,Sz)T  ■  A tDG(l>n+1  =  £>uq-  (S?,STJ,S2)T-cbc  (53) 

Uqn+1  =  fiq  _  (S? ,  S’7,  S*)T  •  AtG4>n+1  (54) 

where  the  T  superscript  denotes  the  transpose.  The  vector  of  volume  fluxes  uq  is  the  contravariant 
component  of  u  in  the  coordinate  system  defined  by  S9  and  Sq. 


3.2.2  Compressible  flow  solver 


Solutions  of  the  compressible  Navier-Stokes  equations  are  obtained  using  Cobaltgo*  an  unstruc¬ 
tured  finite- volume  method  developed  for  solution  of  the  compressible  Navier-Stokes  equations 
and  described  in  Strang  et  al.  (1999).  The  governing  equations  for  a  compressible  flow  written  in 


integral  form  are  expressed  as, 


IJ  QdV  +  JJ  ( fi  +  gj  +  hk)  ■  MS  =  JJ  (n  +  sj  +  tk)  ■  MS ,  (55) 


where, 


’  p " 

pu 

pv 

pw 

pu 

pu 2 

puv 

puw 

pv 

f  = 

puv 

9  - 

pv2  +p 

h  = 

■  pvw 

pw 

puw 

pvw 

pw 2  +p 

-Pe- 

u(pe  +p)_ 

y(pe  +p)_ 

w(pe  +p)_ 

'  o ' 

‘  o  * 

'  0  ' 

Txx 

rxy 

Txz 

Txy 

s  — 

Tyy 

t  = 

Tyz 

TXZ 

Tyz 

Tzz 

a 

_  b  _ 

C 

(56) 


(57) 


In  the  above,  a  =  utxx  +  vrxy  +  wtxz  +  kTx,  b  =  urxy  +  VTyy  +  wryz  +  kTy,  and  c  =  utxz  + 
vryz  +  wtzz  +  kTz.  A  fluid  element  volume  over  which  the  equations  are  enforced  is  denoted  V; 
the  bounding  surface  as  S  with  outward-pointing  unit  normal  n.  The  Cartesian  unit  vectors  are  i, 
j,  and  k.  The  fluid  density  is  denoted  p,  p  is  the  pressure;  u,  v,  and  w  are  the  velocity  components; 
e  is  the  specific  energy  per  unit  volume;  T  is  the  temperature;  k  is  the  thermal  conductivity;  and 
txx,  Tyy,  tzz,  rxy,  Txz ,  and  ryz  are  the  viscous  stress  tensor  components.  The  ideal  gas  law  closes 
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the  system  of  equations  and  the  entire  equation  set  is  non-dimensionalized  by  freestream  density 
and  speed  of  sound. 

Integrating  the  equations  around  finite  volumes  in  the  domain  leads  to  the  semi-discrete  form 
for  the  system, 

Vi—  +  (fMi  +  9M 3  +  hM~k)  ■  hMsM  =  i,  +  sMj  +  tM ’  (58) 

dt  m= 1  m= 1 

where  the  subscripted  i  and  superscripted  M  denote  quantities  for  the  ith  cell  and  the  Mth  face 
of  cell  % ,  respectively,  and  N\  is  the  number  of  faces  bounding  cell  i.  The  equations  above  can 
be  solved  on  arbitrary  cell  types  in  Cobaltfto  (e.g,  hexahedrals,  prisms,  tetrahdrons).  The  spatial 
operator  uses  the  exact  Reimann  solver  of  Gottlieb  and  Groth  (1988),  least  squares  gradient  calcu-. 
lations  using  QR  factorization  to  provide  second  order  accuracy  in  space,  and  TVD  flux  limiters  to 
limit  extremes  at  cell  faces.  A  point  implicit  method  using  analytic  first-order  inviscid  and  viscous 
Jacobians  is  used  for  advancement  of  the  discretized  system.  For  time-accurate  computations,  a 
Newton  sub-iteration  scheme  is  employed,  the  method  is  second  order  accurate  in  time.  Additional 
details  on  the  method  can  be  found  in  Strang  et  al  (1999). 

For  simulation  of  turbulent  flows,  an  averaged  set  of  the  above  equations  are  considered,  this 
process  leading  to  turbulent  stresses  that  must  be  modeled  in  terms  of  other  variables.  A  Boussi- 
nesq  approximation  is  invoked  in  the  momentum  equations  and  the  turbulent  eddy  viscosity,  fit,  is 
used  to  relate  the  stresses  to  the  strain  rate.  The  turbulent  heat  flux  is  also  modeled  using  a  gradient- 
transport  hypothesis,  requiring  specification  of  a  turbulent  thermal  conductivity  kt.  Reynolds  anal¬ 
ogy  is  applied  and  the  turbulent  heat  flux  is  modeled  using  a  constant  turbulent  Prandtl  number, 
Prt  =  {cpfit)/kt,  of  0.9  Using  turbulent  eddy  viscosity  and  turbulent  conductivity,  the  variable  fi 
is  replaced  by  (fi  +  (it)  and  k  is  replaced  by  (k  +  kt)  in  the  governing  equations. 

For  parallel  performance,  Cobalt60  uses  the  domain  decomposition  library  ParMETIS  (Karypis  et 
al  1997)  to  provide  nearly  perfect  load  balancing  with  a  minimal  surface  interface  between  zones. 
Communication  between  processors  is  achieved  using  Message  Passing  Interface  (MPI),  with  par¬ 
allel  efficiencies  above  95%  on  as  many  as  1024  processors  (see  also  Grismer  et  al  1998). 

3.3  Flow  configuration  and  simulation  design 

In  the  first  part  of  the  study,  two-dimensional  computations  of  the  flow  around  canonical  forebody 
cross-sections  were  performed:  a  rounded-comer  square  and  a  2:1  ellipse.  Views  of  both  cross 
sections  for  the  two-dimensional  computations  are  shown  in  Figure  2.  For  the  rounded-comer 
square,  the  radius  is  1/4  of  the  width/height  (“diameter”,  D)  of  the  forebody,  similar  to  the  cross- 
sections  of  the  X-29  and  T-38.  With  reference  to  Figure  2,  the  flow  is  at  angle  of  attack  a  of  10° 
to  the  horizontal  (or)  axis.  This  is  the  angle  of  attack  in  the  idealized  two-dimensional  problem, 
not  the  angle  of  attack  of  the  airplane.  In  simplest  terms,  the  computation  emulates  that  of  an 
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forebody  at  90°  angle  of  attack  (an  exactly  flat  spin),  and  rotating.  The  angle  of  attack  here  is  then 
the  arctangent  of  the  ratio  VtD/V  where  fi  is  the  spin  rate,  D  the  distance  from  the  nose  to  center 
of  rotation  of  the  spin,  and  V  its  sink  velocity. 

For  the  three-dimensional  computations  the  numerical  predictions  of  the  flow  over  the  rounded- 
comer  square  are  compared  to  the  experimental  measurements  from  Polhamus  et  al  (1959).  These 
investigators  measured  the  forces  and  pressure  distributions  on  a  variety  of  forebody  cross-sections 
over  a  range  of  Reynolds  numbers  and  angles  of  attack,  showing  that  there  are  important  Reynolds 
number  effects  in  the  flow  around  forebodies  at  angle  of  attack.  In  the  sub-critical  regime  (Reynolds 
numbers  below  about  5  x  105),  boundary  layer  separation  along  the  top  surface  of  the  rounded- 
comer  square  (upper-most  horizontal  surface  in  Figure  2a)  occurs  near  the  upper-front  comer  of 
the  forebody,  while  for  the  super-critical  flows  the  boundary  layer  remains  attached  along  the  lower 
and  upper  surfaces.  For  the  ellipse  the  change  in  separation  type  is  less  pronounced  given  the  con¬ 
figuration,  with  turbulent  separation  slightly  aft  of  the  upper  shoulder  and  laminar  boundary  layer 
separation  upstream  of  the  upper  shoulder  of  the  model.  The  changes  in  boundary  layer  separation 
characteristics  have  significant  effects  on  the  streamwise  and  lateral  forces  with  Reynolds  number 
(the  lateral,  or  side,  force  acts  along  the  y  axis  in  Figure  2).  A  reversal  of  the  lateral  force  was 
measured  in  the  experiments,  i.e.,  negative  for  sub-critical  flows  and  positive  in  the  super-critical 
regime  (for  reference,  the  force  is  positive  for  an  airfoil).  Relevant  to  spin,  the  negative  side  force 
in  the  sub-critical  regime  is  spin-propelling,  while  at  the  higher  Reynolds  numbers  the  positive  side 
force  is  spin-damping. 

The  effects  measured  by  Polhamus  et  al.  (1959)  are  analogous  to  that  of  the  drag  crisis  oc¬ 
curring  around  cylinders  and  spheres  and  are  linked  to  boundary  layer  transition  and  the  nature  of 
the  flow  separation.  While  predicting  details  of  the  transition  process  is  beyond  the  scope  of  the 
methods  used  in  the  current  simulations,  it  is  possible  to  construct  well-defined  computations  to 
investigate  the  effect  of  the  type  of  boundary  layer  separation  on  the  flow.  In  particular,  simula¬ 
tions  were  performed  in  the  sub-  and  super-critical  regimes  in  which  the  type  of  boundary  layer 
separation  was  controlled  via  the  initial  and  boundary  conditions  on  the  eddy  viscosity. 

Following,  Travin  et  al  (2001),  a  ‘tripless’  approach  is  employed  for  sub-critical  flows  in 
which  the  inflow  eddy  viscosity  is  zero.  Non-zero  values  are  seeded  into  the  wake  and  the  reversed 
flow  that  is  established  in  the  recirculating  region  is  sufficient  to  sustain  the  turbulence  model 
downstream  of  separation.  Boundary  layer  separation  in  this  case  is  laminar,  with  the  model  active 
following  separation.  For  the  super-critical  flows  the  inflow  eddy  viscosity  is  set  to  a  small  value 
(37/),  sufficient  to  ignite  the  turbulence  model  near  solid  surfaces  as  the  fluid  enters  the  boundary 
layers.  The  subsequent  separation  is  then  of  a  turbulent  boundary  layer. 

Computations  of  the  two-dimensional  configurations  presented  below  were  performed  using 
structured  grids,  most  of  the  calculations  of  the  three-dimensional  flow  were  also  performed  on 
structured  meshes  Structured  grids  were  generated  using  the  control  technique  of  Hsu  and  Lee 
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(1991).  Using  this  technique,  it  was  possible  to  control  mesh  spacing  to  the  first  point  nearest  the 
boundary  (within  one  wall  unit  near  solid  surfaces),  exert  control  over  grid  spacing  tangential  to 
the  boundary,  and  also  to  maintain  orthogonality  of  the  mesh  at  all  boundaries.  Important  to  the 
application  of  DES  as  a  tool  for  predicting  the  flow  around  complex  configurations  is  application 
of  the  method  on  unstructured  grids.  For  the  three-dimensional  computations,  DES  predictions  ob¬ 
tained  on  an  unstructured  mesh  were  also  evaluated  against  experimental  measurements  and  results 
from  the  structured  grids.  The  unstructured  grid  was  generated  using  Gridgen  (Steinbrenner  et  al 
2000),  with  prisms  in  the  boundary  layer  and  near-isotropic  tetrahedra  away  from  solid  surfaces. 

Calculations  of  the  three-dimensional  flow  around  the  forebody  cross  section  were  performed 
using  both  structured  and  unstructured  grids.  The  unstructured  grids  are  comprised  of  a  combina¬ 
tion  of  tetrahedra  and  prisms,  while  the  structured  grids  are  comprised  of  hexahedra.  Prisms  are 
used  in  the  boundary  layer  in  order  to  reduce  the  number  of  cells  as  well  as  to  improve  the  effi¬ 
ciency  of  the  boundary  layer  solution.  Boundary-layer  grids  comprised  of  tetrahedra  often  possess 
high  aspect  ratios  and  can  be  strongly  non-orthogonal.  This  presents  problems  in  calculation  of  the 
divergence  of  the  gradient.  Prisms  are  more  orthogonal  and  place  less  burden  on  the  solver. 

Finally,  as  also  discussed  above,  one  of  the  main  objectives  of  the  work  are  to  understand 
forebody  spin  characteristics  as  they  relate  to  side-force  reversal  and  to  assess/advance  DES  as  a 
viable  method  for  prediction  of  unsteady  flows  at  high  Reynolds  numbers.  The  stability  of  DES 
results  with  changes  in  grid  spacing  is  investigated,  as  well  as  other  factors  such  as  the  dimension  of 
the  (statistically  homogeneous)  span  wise  coordinate.  DES  results  are  also  compared  to  predictions 
obtained  from  the  unsteady  Reynolds-averaged  Navier-Stokes  (URANS)  equations  (of  both  the 
two-  and  three-dimensional  equations)  and  to  simulations  performed  without  an,  explicit  turbulence 
model.  The  runs  without  a  turbulence  model  are  denoted  as  MILES  (Monotone  Integrated  Large 
Eddy  Simulation)  to  provide  a  link  with  relevant  literature,  although  no  detailed  investigations  were 
undertaken  to  evaluate  the  numerical  dissipation  in  the  current  calculations  and  its  role  as  an  SGS 
model,  and  the  numerical  schemes  are  not  monotone  in  a  strict  sense  (see  Fureby  and  Grinstein 
1999  for  a  discussion  of  the  MILES  approach). 

4  Results 

4.1  Two-dimensional  computations  -  URANS 

The  initial  phase  of  the  work  was  undertaken  to  develop  and  build  confidence  in  the  approach  used 
to  specify  and  control  the  type  of  boundary  layer  transition,  i.e.,  either  fully  turbulent  boundary 
layer  separation  or  the  tripless  approach  used  for  cases  with  laminar  boundary  layer  separation. 
Unsteady  RANS  comprised  the  primary  simulation  technique  for  this  stage  of  the  work,  the  S- 
A  model  was  employed  to  compute  the  eddy  viscosity  used  to  close  the  Reynolds  stresses.  The 
numerical  solutions  were  of  the  unsteady  incompressible  flow  using  the  fractional  step  method 
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outlined  in  §3.2.1.  This  portion  of  the  study  also  enabled  implementation  and  testing  of  the  grid 
control  strategy  developed  by  Hsu  and  Lee  (1991)  for  generation  of  structured  meshes  about  the 
forebodies.  A  grid  generation  procedure  was  developed  and  tested  for  two  canonical  forebody 
cross-sections  described  above:  a  rounded-comer  square  and  a  2:1  ellipse.  Shown  in  Figures  2-4 
are  views  of  each  geometry.  In  Figure  3  and  Figure  4  show  the  grid  and  entire  computational 
domain  (frame  a  of  each  figure),  a  view  in  the  vicinity  of  the  forebody  is  also  shown  (frame  b  of 
each  figure). 

Important  in  application  of  the  S-A  model  for  URANS,  as  well  as  in  the  DES  predictions 
presented  in  the  next  section,  is  resolution  of  the  wall  layer,  i.e.,  the  first  grid  point  nearest  a 
solid  surface  is  located  approximately  one  viscous  unit  from  the  surface  with  geometric  stretching 
applied  in  the  boundary  layer.  The  stretching  rate  applied  to  the  boundary  layer  grids  were  in  the 
range  1.2  -  1.3,  with  the  higher  value  representing  a  very  aggressive  stretching.  The  position  of 
the  first  wall-normal  grid  point  is  estimated  via, 


where  rw  is  the  wall  stress  and  related  to  the  skin  friction  via  Cf  =  2 Tw/(pU^).  The  above 
expression  can  be  rewritten  in  terms  of  the  skin  friction  coefficient  Cf , 

(60) 

where  D  is  the  body  width/diameter  and  Re&  =  UooD/u.  An  estimate  of  y\  corresponding  to  a 
viscous  spacing  approximately  one  wall  unit  can  be  obtained  using  a  representative  value  of  C /, 
e.g,  in  the  range  Cf  «  0.002  -  0.004,  the  higher  value  yielding  a  more  conservative  estimate  of  the 
near-wall  spacing  to  the  first  point.  For  y+  =  1,  the  formula  (60)  yields  yi/D  in  the  range  10-5  for 
the  Reynolds  numbers  considered  in  this  work,  Rep  on  the  order  105  -  106. 

For  the  rounded-comer  square,  Figure  3  shows  the  axial  distance  of  the  inlet  from  the  front 
of  the  forebody  is  3.5D,  the  downstream  extent  7.5 D,  with  the  lateral  dimension  extending  3.5D 
at  the  maximum  distance  from  the  forebody.  The  domain  for  the  2:1  ellipse  is  similar,  i.e.,  using 
the  long  axis  as  the  characteristic  length,  the  streamwise  coordinate  ranges  —4 D  <  x  <  8 D, 
the  lateral  dimension  -4D  <  y  <  4 D.  Figure  3b  and  Figure  4b  show  that  near  the  forebody 
surface,  the  wall-normal  grids  are  clustered,  in  the  surface-tangent  dimension  the  grids  are  shaded 
towards  the  wake  with  a  higher  density  of  cells  in  the  aft  region.  Control  of  these  features  in 
the  grid  are  relatively  straightforward  to  achieve  using  the  approach  developed  by  Hsu  and  Lee 
(1991),  a  further  advantage  being  the  relatively  low  cell  skewness.  Grid  refinement  was  carried  out 
using  simulations  in  which  the  surface-normal  coordinate  was  refined  from  100  to  150  points,  the 
surface-tangential  dimension  was  refined  from  150  to  200  points. 

Calculations  of  the  two-dimensional  configurations  were  performed  for  Reynolds  numbers  of 
5  x  104  and  8  x  105.  The  lower  Reynolds  number  was  employed  for  the  initial  investigations  of 
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the  tripless  characteristics  of  the  model,  simulations  at  this  Reynolds  number  lowering  the  com¬ 
putational  burden  (e.g.,  allowing  slightly  coarser  meshes)  and  also  alleviating  difficulties  initially 
encountered  with  the  tripless  approach.  Eddy  viscosity  levels  were  observed  to  “creep”  upstream 
of  separation  due  to  numerical  errors,  this  aspect  necessitating  the  modification  to  the  production 
term  using  fv3  as  described  in  §3.1.  An  example  of  the  influence  of  boundary  layer  separation  on 
the  nature  of  the  solution  along  the  top  surface  of  the  rounded-comer  square  is  shown  in  Figure  5. 
Both  frames  of  the  figure  show  instantaneous  velocity  vectors  colored  by  the  eddy  viscosity  ratio 
vtjv .  The  computations  corresponding  to  each  frame  were  performed  Re&  =  5  x  104,  in  Figure  5a, 
the  boundary  layers  on  the  forebody  surface  are  fully  turbulent,  this  feature  accomplished  by  seed¬ 
ing  the  inlet  plane  with  a  small  level  of  eddy  viscosity,  vt/v  —  3.  This  level  of  eddy  viscosity  is 
sufficient  to  activate  the  turbulence  model  in  the  boundary  layers,  a  consequence  and  illustrated 
in  Figure  5a  is  that  the  boundary  layers  possess  sufficient  momentum  to  withstand  the  adverse 
pressure  encountered  as  the  flow  develops  along  the  top  surface  of  the  forebody.  Consequently, 
boundary  layer  separation  is  delayed,  occurring  near  the  upper  rear  comer  as  shown  in  the  figure. 
Figure  5b  shows  instantaneous  velocity  vectors  from  the  tripless  solution.  The  tripless  solution  is 
created  by  seeding  the  initial  flow  with  eddy  viscosity  and  with  zero  eddy  viscosity  prescribed  at 
the  inlet  to  the  computational  domain.  The  reverse  flow  region  sustains  the  model,  Figure  5b  shows 
that  boundary  layer  separation  occurs  near  the  upper  front  comer  of  the  forebody.  In  this  case,  the 
boundary  layers  developing  along  the  forebody  are  effectively  laminar  and  do  not  possess  sufficient 
momentum  to  withstand  the  adverse  pressure  gradient  that  develops  as  the  flow  progresses  around 
the  upper  front  comer  and  along  the  top  surface.  The  larger  region  of  non-zero  eddy  viscosity  (red 
vectors  in  the  figure)  is  due  to  the  separation,  the  recirculating  region  that  has  developed  along  the 
upper  surface  sweeping  eddy  viscosity  into  the  regions  above  the  forebody.  Importantly  for  this 
case,  the  turbulence  model  remains  dormant  in  the  region  upstream  of  separation. 

Variations  in  the  pressure  coefficient,  skin  friction,  and  “turbulence  index”  around  the  perimeter 
of  the  2:1  ellipse  are  shown  in  Figure  6.  The  distributions  are  shown  for  a  tripless  solution  at 
ReD  =  5  x  104  in  Figure  6a  and  a  case  with  fully  turbulent  boundary  layers  at  ReD  =  8  x  105 
in  Figure  6b.  The  horizontal  coordinate  defines  the  angle  around  the  perimeter  and  is  measured 
positive  in  the  clockwise  direction.  The  coordinates  6  =  0  and  360  degrees  correspond  to  the 
windward  symmetry  point,  the  upper  and  lower  “shoulders”  of  the  ellipse  coinciding  with  6  —  90 
and  270  degrees  (where  the  potential  flow  at  0  degrees  angle  of  attack  would  possess  the  minimum 
pressure).  The  pressure  coefficient  distribution  in  both  Figure  6a  and  Figure  6b  possess  common 
features  on  the  front  of  the  body,  the  stagnation  pressure  is  predicted  near  9  =  0  and  there  is 
a  relatively  rapid  decrease  as  the  flow  accelerates  around  the  upper  and  lower  shoulders  of  the 
ellipse.  The  back  pressure  for  the  case  with  laminar  separation  in  Figure  6a  shows  larger  suction, 
an  effect  that  will  lead  to  a  larger  axial  force.  The  skin  friction  distribution  for  both  cases  shows  the 
zero  crossing  marking  the  location  of  flow  separation.  For  the  2:1  ellipse,  separation  for  the  two 
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cases  occurs  at  similar  angles,  only  slightly  delayed  for  the  fully  turbulent  case  around  the  lower 
corner,  occurring  at  nearly  the  same  angle  around  the  upper  comer.  For  both  the  tripless  and  fully 
turbulent  solutions,  Cf  increases  as  the  flow  accelerates  towards  the  upper  and  lower  shoulders  and 
then  is  nearly  zero  in  the  separated  region. 

Also  shown  in  the  figure  is  the  “turbulence  index”,  which  defines  the  position  at  which  the 
turbulence  model  becomes  active.  For  the  S-A  model,  the  turbulence  index  is  defined  as, 


dv 
dn  5 


(61) 


where  the  coordinate  n  defines  the  surface  normal.  The  index  is  zero  in  a  laminar  region,  followed 
by  a  sharp  increase  in  the  “transition  region”,  i.e.,  the  region  over  which  the  model  becomes  active, 
and  then  remains  around  unity  or  takes  on  higher  values  in  the  fully  turbulent  regions.  Figure  6a 
shows  that  for  the  tripless  case  the  abrupt  increase  in  the  turbulence  index  coincides  with  the 
zero  crossing  in  the  skin  friction,  i.e.,  with  boundary  layer  separation,  in  turn  showing  the  tripless 
character  of  the  solution  in  that  the  model  remains  dormant  in  the  attached  boundary  layer  upstream 
of  separation.  Figure  6b  shows  that  for  the  fully  turbulent  boundary  layers  the  turbulence  index  is 
essentially  equal  to  unity  prior  to  separation,  abruptly  increasing  as  the  flow  separates  (the  sharp 
increase  due  to  small/zero  values  in  uT). 

An  example  of  the  unsteady  character  of  the  URANS  solution  is  illustrated  for  the  rounded- 
comer  square  in  Figure  7.  Shown  are  vorticity  contours  for  several  phases  of  the  motion  through 
roughly  one  shedding  cycle.  The  calculations  are  from  a  fully  turbulent  case  at  ReD  =  8  x 
105  for  which  the  turbulent  boundary  layers  separate  near  the  upper  and  lower  surfaces  of  the 
model,  an  effect  observed  in  the  contours.  The  vorticity  contours  show  essentially  a  periodic 
shedding,  with  well  organized  structures  developing  in  the  wake,  the  figure  also  showing  that 
the  asymmetry  that  develops  in  the  wake  due  to  the  10°  angle  of  attack  of  the  freestream.  The 
distributions  of  the  pressure  coefficient,  skin  friction,  and  turbulence  index  corresponding  to  the 
frames  shown  in  Figure  7  are  shown  in  Figure  8.  As  observed  earlier  for  the  ellipse,  the  turbulence 
index  is  approximately  equal  to  one  in  the  attached  boundary  layers  upstream  of  separation,  sharply 
rising  following  flow  detachment.  Zero  crossings  in  the  skin  friction  again  identify  the  separation 
location.  For  the  fully  turbulent  solution,  Figure  7  shows  that  separation  occurs  near  the  rear  of 
the  upper  and  lower  horizontal  surfaces,  corresponding  to  angles  8  ~  130°  (top  surface)  and  9  « 
230°  (bottom  surface).  For  all  of  the  phases  of  the  shedding,  the  pressure  coefficient  distribution 
shows  that  there  is  a  very  large  acceleration  as  the  flow  negotiates  the  upper  front  comer,  the 
deep  decrease  in  Cv  occurring  because  the  boundary  layer  remains  attached  through  this  region. 
A  smaller  minima  in  Cp  occurs  near  the  lower  front  comer  ( 8  «  315°)  due  to  boundary  layer 
acceleration.  The  skin  friction  distribution  shows  that  there  is  no  boundary  layer  separation  along 
the  lower  horizontal  surface.  A  clearly  observed  effect  of  vortex  shedding  is  the  change  in  the 
pressure  distribution  along  the  rear  surface  of  the  forebody,  for  130  <  8  <  225  degrees.  Figure  8a, 
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Case 

Model 

Grid  Size 

u. 

x-y  domain 

i 

DES 

100  x  149  x  151 

3  D 

baseline 

2 

DES 

100  x  149  x  301 

6D 

baseline 

3 

DES 

150  x  200  x  151 

3  D 

baseline 

4 

DES 

120  x  149  x  151 

3  D 

padded 

5 

DES 

unstructured,  3.55  x  106  cells 

3  D 

padded 

6 

MILES 

120  x  149  x  151 

3D 

padded 

7 

URANS 

200  x  400 

- 

baseline 

8 

URANS 

120  x  149  x  151 

3D 

padded 

Table  1:  Simulation  parameters.  Grid  size  reported  as  surface-normal  x  surface-tangential  x 
spanwise;  “baseline”  is  the  smaller  x—y  domain,  “padded”  the  larger  x-y  domain.  URANS  calcu¬ 
lations  performed  using  the  Spalart-Allmaras  one-equation  model. 


for  example  shows  a  relatively  large  suction  along  the  rear  vertical  surface,  the  corresponding  flow 
visualization  in  Figure  7a  shows  a  vortex  rolling  up  directly  behind  the  rear  vertical  surface. 

4.2  Three-dimensional  computations  -  DES  and  URANS 

In  the  first  part  of  the  study  the  tripless  approach  to  producing  solutions  undergoing  laminar  bound¬ 
ary  layer  separation  and  with  a  turbulent  wake  were  developed  and  tested.  Overall,  the  behavior 
of  the  approach  is  well-defined,  not  sensitive  to  initial  or  boundary  conditions,  etc.  In  the  second 
phase  of  the  work,  three-dimensional  computations  of  the  flow  around  the  forebody  were  per¬ 
formed,  a  schematic  of  the  geometry  and  indication  of  the  angle  of  the  incident  flow  is  shown  in 
Figure  9  for  convenience. 

The  parameters  of  the  three-dimensional  calculations  are  summarized  in  Table  1.  Shown  is  the 
case  number,  model,  grid  size,  spanwise  period,  and  reference  to  the  x—y  domain  size.  Assessment 
of  unstructured  approaches  to  eddy-resolving  calculations  is  enabled  via  Case  5,  the  ability  to 
exert  greater  control  on  cell  distribution  compared  to  the  structured  grids  permitted  generation  of 
an  unstructured  mesh  having  2.5  x  106  cells  (of  a  total  cell  size  of  3.55  x  106  cells)  within  two 
diameters  of  the  model  surface. 

The  majority  of  simulations  were  performed  on  domains  in  which  the  spanwise  dimension  was 
three  times  the  diameter,  i.e.,  Lz  =  3 D.  The  influence  of  the  spanwise  period  was  investigated  in 
Case  2  via  computations  performed  on  a  domain  with  a  doubling  of  the  spanwise  period  to  Lz  =  6. 
Domain-size  influences  were  also  investigated  in  calculations  with  two  domains  having  different 
outer-boundary  placement.  The  smaller  domain  (referred  to  as  “baseline”  in  the  table)  extended 
eight  diameters  downstream  of  the  body  and  with  a  lateral  extent  also  of  8 D,  a  projection  of  the 
grid  onto  the  x-y  plane  for  the  baseline  configuration  is  shown  in  Figure  10. 

In  calculations  performed  on  the  larger  domains  (referred  to  as  “padded”  in  the  table),  the 
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streamwise  extent  to  the  outflow  boundary  downstream  of  the  body  was  at  approximately  20 D, 
with  the  lateral  dimension  also  approximately  20D  from  the  model  surface.  As  shown  in  the  next 
section,  there  was  a  measurable  effect  of  the  baseline  domain  on  the  solutions,  resulting  in  over¬ 
predictions  of  the  stagnation  pressure  and  axial  force.  Boundary  conditions  on  the  model  surface 
were  no-slip  for  the  velocity  components  and  turbulent  viscosity.  The  normal  momentum  equation 
is  used  at  solid  walls  to  estimate  the  variation  of  pressure  normal  to  the  surface,  while  a  one-sided, 
least-squares  gradient  method  is  used  to  estimate  the  variation  of  pressure  tangential  to  the  wall. 

Grid  resolution  effects  were  investigated  by  refining  the  mesh  in  the  x-y  plane  by  a  factor 
of  two.  For  the  coarse  grid  with  the  shorter  spanwise  dimension  (Lz  =  3D),  the  structured- 
grid  calculation  was  performed  using  approximately  2.2  x  106  points,  the  finer  mesh  calculation 
possessing  about  4.5  x  106  grid  points.  Because  A  in  (12)  near  the  surface  was  set  by  the  spanwise 
spacing,  the  thickness  of  the  “RANS  region”,  i.e.,  the  dimension  from  the  model  surface  to  the 
interface  at  which  d  is  set  by  the  grid  spacing  for  all  the  DES  runs  was  0.013D.  The  dimensionless 
timestep,  At/{D/Uoo)  (Doo  is  the  freestream  speed),  was  0.01,  a  conservative  value  chosen  based 
on  preliminary  calculations  and  previous  time-accurate  computations  of  unsteady  flows  by  the 
current  investigators  and  other  researchers  (e.g.,  see  Travin  et  al  2001).  With  At/(D/Uoo)  =  0.01, 
there  are  approximately  350  timesteps  per  main  shedding  cycle. 

From  a  given  set  of  initial  and  boundary  conditions  for  a  particular  flow  type  (tripless  or  fully- 
turbulent),  the  governing  equations  were  time  advanced  through  a  transient  as  the  flow  evolved 
to  its  equilibrium  condition.  This  transient,  typically  less  than  20D/f7oo»  was  discarded  and  the 
simulations  continued  for  an  additional  period  of  0(100D/C/oo)*  This  period  was  sufficient  for 
adequate  convergence  of  averaged  quantities  and  capture  of  the  long  timescales  in  the  flow.  Initial 
perturbations  supplied  to  the  flow  triggered  the  breakdown  in  the  wake  and  development  of  three- 
dimensional  structures. 

A  snapshot  of  the  instantaneous  velocity  vectors  in  a  plane  near  the  upper  rear  comer  ( 9  in  the 
range  135°)  for  a  turbulent  separation  run  (Case  4)  is  shown  in  Figure  11.  The  interface  beyond 
which  d  is  set  by  the  grid  has  been  drawn.  The  figure  shows  a  smooth  transition  between  the 
“RANS  region”  and  “LES  region”  of  the  solution  with  essentially  all  of  the  boundary  layer  within 
the  RANS  region.  In  addition,  the  outer  part  of  the  boundary  layer  has  very  little  sensitivity  to  the 
destruction  term  of  the  S-A  model,  the  only  one  that  changes  between  RANS  and  LES  modes. 

The  different  structure  of  the  sub-  and  super-critical  flows  is  illustrated  in  Figure  12.  For  the 
sub-critical  flow  at  Re  =  105  the  attached  boundary  layers  are  laminar  and  cannot  sustain  the 
development  of  the  adverse  pressure  gradient  in  the  upper  front  comer  with  separation  occurring 
around  9  «  45°.  In  the  “tripless”  mode,  the  turbulence  model  is  dormant  upstream  of  separation, 
with  flow  reversal  sustaining  the  model  downstream  of  separation.  Contrasted  with  Figure  12b, 
the  tripless  case  exhibits  massive  separation  along  the  upper  surface,  the  fully  turbulent  solution 
in  Figure  12b  shows  the  influence  of  the  fully  turbulent  boundary  layers  along  the  forebody,  in 
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particular  the  attached  flow  along  the  top  surface  with  separation  towards  the  rear.  Both  flows 
show  that  the  wake  is  chaotic  and  three-dimensional. 

Shown  in  Figure  13  are  a  representative  sample  of  pathlines  in  the  vicinity  of  the  upper  front 
comer  of  the  forebody  for  the  flow  at  Re  =  105  (Case  4  in  Table  1).  The  pathlines  in  the  figure  are 
colored  by  the  value  of  the  viscosity  ratio  vt/v.  They  cross  in  this  view,  because  the  flow  is  three- 
dimensional.  Upstream  of  separation  the  eddy  viscosity  is  zero  (as  indicated  by  the  blue  color  of 
the  pathline).  In  the  separated  region  the  reversing  flow  sweeps  turbulent  fluid  from  downstream 
into  contact  with  the  separating  flow.  Important  to  note  is  that  there  is  not  a  “transition  creep”,  i.e., 
a  numerical  diffusion  of  non-zero  eddy  viscosity  into  regions  upstream  of  separation. 

In  “natural”  applications  of  DES,  the  detached  regions  of  the  flow  are  computed  using  LES, 
in  this  case  with  a  one-equation  model  for  the  subgrid-scale  eddy  viscosity.  An  advantage  of  LES 
is  that  mesh  refinement  resolves  more  flow  features,  in  turn  lessening  modeling  errors  and  driving 
the  solution  towards  the  DNS  limit.  The  effect  of  mesh  refinement  was  investigated  by  doubling 
the  x—y  grid  in  Case  3  as  compared  to  Case  1  (c.f.,  Table  1).  Shown  in  Figure  14a  and  Figure  14b 
are  instantaneous  vorticity  contours  from  Case  1  and  Case  3,  respectively,  for  a  simulation  with 
turbulent  separation  at  Re  =  8  x  105.  Cuts  of  the  vorticity  field  from  three  spanwise  planes  are 
shown  for  each  case  and  provide  an  example  of  the  strong  spanwise  variation  in  the  DES  solution. 
As  also  the  case  for  classical  LES,  Figure  14  shows  that  the  effect  of  the  mesh  refinement  is  to 
resolve  smaller- scale  eddies  in  Case  3.  This  feature  in  the  DES  was  also  illustrated  in  the  circular 
cylinder  calculations  of  Travin  et  al  (2001). 

Also  shown  in  Figure  14  are  vorticity  contours  in  the  forebody  wake  from  solutions  obtained 
using  the  unstructured  grid  (Case  5)  and  a  tripless  solution  with  simulation  parameters  correspond¬ 
ing  to  Case  3.  The  unstructured-grid  prediction  shows  roughly  a  comparable  range  of  eddies  re¬ 
solved  in  the  wake  as  the  fine-mesh  result  obtained  for  Case  3.  That  the  unstructured  result  yields 
similar  structure  as  Case  3  is  significant  since  it  demonstrates  that  unstructured  meshes  provide 
an  efficient  means  for  performing  DES,  the  ability  to  to  exert  a  more  local  control  on  the  mesh  in 
the  focus  region  in  the  wake  of  the  body  an  important  aspect  for  turbulence-resolving  simulations. 
Vorticity  contours  from  the  tripless  solution  in  Figure  14d  show  a  substantially  different  structure 
than  achieved  in  the  fully  turbulent  solutions.  Boundary  layer  separation  is  apparent  near  the  upper 
front  comer  of  the  forebody,  with  a  wider  wake  developing  downstream.  A  zoomed  view  of  the 
vorticity  contours  is  shown  in  Figure  15  shows  instantaneous  vorticity  contours  in  the  wake  for 
the  fully  turbulent  solution  from  the  structured-grid  prediction  (Case  3)  and  the  unstructured-grid 
prediction  (Case  5).  Visually,  the  range  of  scales  resolved  in  the  wake  is  comparable  for  the  two 
grids. 

Force  coefficients  Cx  and  Cy  in  the  axial  and  lateral  directions,  respectively,  are  defined  using 
the  freestream  density  and  velocity  and  frontal  area  of  the  forebody.  Figure  16  shows  the  force 
coefficient  time  histories  for  the  two-dimensional  URANS  solution  at  Re  =  8  x  105.  A  fraction  of 
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the  time  history  is  shown  in  the  figure,  i.e.,  following  the  initial  transient.  The  figure  shows  the  2D 
URANS  solution  is  temporally  periodic,  with  large  swings  in  the  side  force  coefficient  compared 
to  the  axial  value  as  the  flow  undergoes  a  shedding  cycle. 

A  representative  force  coefficient  history  for  a  DES  run  is  shown  in  Figure  16,  for  a  turbulent 
separation  run  at  Re  =  8  x  105  (Case  2).  Similar  to  Figure  16,  a  transient  of  roughly  20  non- 
dimensional  time  units  has  been  excluded  from  the  figure  (note  also  the  longer  time  integration  for 
the  DES).  Unlike  the  2D  URANS,  a  strong  modulation  is  apparent  in  the  side  force  coefficient  Cy, 
similar  to  that  observed  in  related  studies  of  cylinders  and  other  bluff  bodies  (e.g.,  see  Travin  et  al. 
2001).  The  side-force  modulation  is  complex  and  seems  to  be  an  intrinsic  feature  of  the  chaotic, 
three-dimensional  flow.  For  the  forebody,  the  modulation  develops  via  the  interaction  of  spanwise 
and  streamwise  vorticity  in  the  near  wake.  DES  calculations  on  domains  in  which  the  spanwise 
coordinate  Lz  was  1.5D  did  not  yield  force  modulation  and  suppressed  three-dimensionality  of 
the  primary  spanwise  structure  (although  the  solution  possessed  streamwise  vorticity).  Predictions 
on  the  domain  with  Lz  -  1.5 D  yielded  large  over-predictions  of  the  axial  force.  Though  not 
shown  here,  force-coefficient  histories  for  all  of  the  three-dimensional  turbulent  separation  cases  - 
including  the  3D  URANS  result  -  exhibited  force  modulation. 

Time-averaged  force  coefficients  for  the  turbulent  separation  cases  are  summarized  in  Table  2. 
The  2D  URANS,  which  produces  a  periodic  shedding  and  cannot  accurately  account  for  the  force 
modulation,  substantially  over-predicts  the  mean  axial  force  coefficient  ( Cx ).  This  feature  is  analo¬ 
gous  to  the  circular  cylinder  where  two-dimensional  URANS  yields  large  drag  (Travin  et  ah  2001). 
For  the  DES,  force  coefficients  from  the  smaller  (baseline)  domains  are  higher  than  the  measured 
values  and  than  those  from  calculations  performed  on  the  larger  domain  (c.f..  Cases  1-3).  A  com¬ 
parison  of  Case  3  against  Cases  1-2  show  a  trend  towards  lower  axial  force  with  grid  refinement 
in  the  x-y  plane.  In  addition,  the  axial  force  slightly  decreases  in  Case  2,  computed  on  the  longer 
spanwise  domain,  compared  to  Case  1.  Nevertheless,  ( Cx )  is  too  high  and,  as  shown  below,  the 
over-prediction  arises  from  the  influence  of  the  computational  domain,  which  effectively  constrains 
the  flow  and  raises  the  stagnation  pressure  coefficient  by  about  0.1  compared  to  results  obtained  on 
the  larger  domains  (Cases  4-5).  DES  predictions  on  the  larger  domain  using  both  structured  and 
unstructured  grids  are  in  quite  good  agreement  with  the  measurements  of  Polhamus  et  al.  (1959). 
Interestingly,  the  3D  URANS  also  yields  forces  in  good  agreement  with  the  measurements  and 
padded-domain  DES  predictions.  On  the  other  hand,  calculations  without  an  explicit  turbulence 
model  (Case  6)  markedly  over-predict  the  axial  force  due  to  the  poor  treatment  of  the  attached 
boundary  layers,  as  described  in  more  detail  below. 

Pressure  coefficients  around  the  body  for  the  fully -turbulent  runs  (Re  =  8  x  105)  are  shown 
in  Figure  17.  The  angle  6  is  measured  counter-clockwise  from  the  aft-symmetry  point  of  the 
forebody.  For  the  flow  at  10°  angle  of  attack,  the  maximum  Cp  occurs  about  15  -  20°  below  the 
fore-symmetry  point  (6  ~  —160°  as  shown  in  Figure  17).  For  Cases  1,  3,  and  7,  the  stagnation 


24 


Case 

Model 

(Cy) 

i 

DES 

0.57 

0.92 

2 

DES 

0.55 

0.98 

3 

DES 

0.51 

0.96 

4 

DES 

0.46 

0.94 

5 

DES 

0.43 

0.83 

6 

MILES 

0.76 

0.62 

7 

2D  URANS 

0.75 

0.88 

8 

3D  URANS 

0.43 

0.94 

- 

expts. 

0.4 

0.9 

Table  2:  Time-average  force  coefficients  from  turbulent  separation  cases  at  Re  —  8  x  105  (time 
averages  denoted  using  {•}).  Experimental  measurements  are  from  Polhamus  et  al  (1959). 


Cp  is  over-predicted,  an  error  introduced  by  the  use  of  the  smaller  x-y  domain.  For  these  cases 
the  over-prediction  in  the  stagnation  pressure  is  (9(0.1),  comparable  to  the  over-prediction  of  the 
mean  axial  force  (c.f.,  Table  2).  Comparing  the  effect  of  the  domain  (Case  1  and  Case  4)  shows 
a  reduction  in  the  axial  force  coefficient  on  the  larger  domain,  closer  to  the  measured  value  of 
about  0.4  as  also  summarized  in  Table  2.  The  2D  URANS,  with  a  finer  x-y  grid  compared  to  the 
three-dimensional  runs  also  shows  deeper  minima  in  Cv  in  the  vicinity  of  each  comer,  providing 
some  insight  into  the  effect  of  grid  resolution  on  the  pressure  field. 

An  effect  of  the  Reynolds  number  reproduced  in  the  fully-turbulent  solutions  in  Figure  17  is 
that  the  boundary  layers  around  the  upper  front  and  lower  front  comers  ( 9  ^  ±135°)  remain  at¬ 
tached,  as  evidenced  by  the  strong  pressure  minima  in  these  regions,  especially  around  9  ~  135°. 
It  is  apparent  that  all  of  the  simulations  predict  attached  boundary  layers  around  the  upper  front 
comer,  with  the  exception  of  the  MILES  run  (Case  6),  i.e.,  the  simulation  performed  without  an 
explicit  turbulence  model  in  which  the  (laminar)  boundary  layer  separates.  Accurate  prediction  of 
boundary  layer  growth  and  separation  in  MILES  requires  the  boundary  layer  grid  be  sufficiently 
fine  to  resolve  the  small  near-wall  turbulent  structures  (as  would  also  be  required  in  whole-domain 
LES  with  an  explicit  SGS  model).  In  practice,  however,  boundary  layer  grids  will  not  be  suffi¬ 
ciently  dense  for  high  Reynolds  numbers  flows  because  of  the  high  computational  cost.  In  the 
present  case,  the  resulting  MILES  boundary  layer  prediction  lacks  turbulent  structures  and  is  es¬ 
sentially  laminar. 

Along  the  rear  vertical  surface  (in  the  vicinity  9  =  0  in  Figure  17)  the  2D  URANS  (Case  7) 
pressure  distribution  is  far  from  the  measured  values,  resulting  in  a  large  stream  wise  force  (c.f.,  Ta¬ 
ble  2)  as  also  observed  in  other  two-dimensional  URANS  predictions  of  bluff-body  flows  (e.g.,  see 
Travin  et  al  2001).  For  all  other  cases  shown  in  the  figure,  including  the  3D  URANS,  predictions 
of  the  rear-surface  pressures  are  reasonable,  close  to  the  measurements  of  Polhamus  et  al  (1959). 
Consequently,  for  the  DES  and  3D  URANS  calculations  on  the  padded  domains  (Cases  4,  5,  and 
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8),  the  overall  pressure  distributions  are  adequate  and  the  axial  and  streamwise  forces  exhibit  rea¬ 
sonable  agreement  with  measurements.  The  accuracy  of  the  3D  URANS  result  is  surprising  since 
flow  visualizations  show  the  solution  lacks  the  streamwise  vorticity  apparent  in  the  DES  predic¬ 
tions  shown  in  Figure  12.  Flow  visualizations  show  that  the  3D  URANS  solution  exhibits  weak, 
but  persistent,  three-dimensionality  in  the  primary  spanwise  vorticity  shed  from  the  forebody.  Be¬ 
cause  the  peak  suction  is  missed  in  the  MILES  run  (Case  6),  the  axial  force  is  too  high,  yielding  a 
similar  (Cx)  to  the  2D  URANS,  though  the  causes  for  the  over-predictions  by  these  two  techniques 
are  not  the  same. 

Mean  side-force  coefficients  summarized  in  Table  2  show  that  the  DES  predictions  of  the  lateral 
force  coefficient  (Cy)  are  in  general  not  far  from  the  measurements  reported  by  Polhamus  et  al 
(1959).  The  lateral  force  prediction  in  the  MILES  case  provides  another  illustration  of  the  error  that 
can  arise  due  to  the  boundary  layer  treatment  in  this  technique.  As  also  noted  in  the  Cv  distribution 
and  axial  force  coefficient,  the  3D  URANS  is  again  accurate  and  apparently  able  to  resolve  enough 
of  the  3D  variation  important  to  accurate  force  predictions. 

Laminar  separation  cases  in  the  DES  were  computed  for  most  of  the  parameter  combinations 
summarized  in  Table  L  The  Reynolds  number  in  the  laminar  separation  runs  was  1  x  105  for  the 
DES  (4  x  105  for  the  2D  URANS  shown  below),  well  below  the  critical  value  found  by  Polhamus  et 
al  (1959)  of  approximately  5  x  105.  In  the  tripless  mode,  the  upstream  eddy  viscosity  was  zero. 
The  wake  was  initially  seeded  with  eddy  viscosity  and  the  reversing  flow  established  behind  the 
forebody  is  sufficient  to  sustain  the  turbulence  model  following  separation. 

A  representative  force  history  from  a  laminar  separation  case  is  shown  in  Figure  18.  The  sim¬ 
ulation  parameters  for  this  case  correspond  to  those  of  Case  4  in  Table  L  As  shown  qualitatively 
via  Figure  12,  below  the  critical  Reynolds  number  the  flow  separates  in  the  vicinity  of  the  upper 
front  comer  of  the  forebody.  The  measurements  of  Polhamus  et  al  (1959)  indicate  that  the  bound¬ 
ary  layer  along  the  lower  surface  of  the  forebody  remains  attached.  The  pressure  distribution  then 
develops  lower  pressures  along  the  lower  forebody  surface  compared  to  the  upper  surface,  which 
has  the  result  of  reversing  the  magnitude  of  the  side  force  as  compared  to  the  values  measured  at 
higher  Reynolds  numbers,  past  the  critical  value. 

The  force  histories  shown  in  Figure  18  show  a  higher  axial  force  than  in  the  fully  turbulent 
runs.  The  mean  axial  force  for  this  case  is  around  0.8,  not  far  from  the  value  reported  in  Table  2 
for  the  MILES  run  at  Re  =  8  x  105,  which  also  experiences  flow  detachment  in  approximately 
the  same  region  near  the  upper  front  comer.  More  importantly,  the  side  force  Cy  in  Figure  18, 
while  chaotic,  is  only  infrequently  negative.  Therefore,  the  mean  side  force  will  not  be  negative 
(the  mean  Cy  is  0.38  for  the  trace  in  Figure  18)  and  the  simulation  does  not  yield  a  reversal  in  the 
magnitude  of  the  side  force.  Changing  the  type  of  separation  produces  a  very  tangible  move  in  the 
correct  direction,  but  still  an  insufficient  one. 

The  pressure  coefficient  distribution  around  the  forebody  for  the  laminar  separation  DES  (Case  4 
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parameters)  is  shown  along  with  the  experimental  measurements  of  Polhamus  et  al  (1959).  For 
comparison,  the  Cp  distribution  from  the  2D  URANS  calculation  at  Re  —  4  x  105  is  also  shown. 
Consistent  with  the  flow  visualization  shown  in  Figure  12,  flow  detachment  around  0  «  135°  re¬ 
sults  in  a  substantially  higher  Cp  compared  to  the  turbulent  separation  case  in  Figure  17.  Both 
the  DES  and  2D  URANS  have  lower  minima,  indicative  of  delayed  boundary  layer  separation  as 
compared  to  the  experiments.  Near  the  lower  front  surface  ( 9  &  —135°),  measurements  show  a 
deeper  minimum  than  predicted  in  the  DES.  The  2D  URANS,  on  the  other  hand,  comes  closer  to 
predicting  the  pressure  minima  along  the  lower  surface.  Along  the  rear  vertical  surface  ( 8  «  0°), 
DES  predictions  of  the  pressure  distribution  are  reasonable.  However,  because  of  the  deeper  min¬ 
ima  in  the  DES  Cv  near  0  «  135°  and  higher  Cv  along  the  lower  flat  surface,  side  force  reversal 
cannot  occur.  Inspection  of  the  instantaneous  fields  shows  that  along  the  lower  flat  surface  (in  the 
vicinity  of  9  «  135°),  a  thin  region  of  reversed  flow  occurs  in  the  DES  (in  the  mean).  This  reversed 
flow  region  contributes  to  an  effectively  altered  geometry  that  prevents  development  of  a  deep  Cp 
minimum  as  apparently  occurs  in  the  experiments.  To  determine  if  the  development  of  the  thin 
region  of  reversed  flow  was  caused  by  numerical  and/or  modeling  errors,  a  Direct  Numerical  Sim¬ 
ulation  of  the  two-dimensional  flow  at  Re  =  1  x  104  was  conducted  using  a  grid  of  1600  x  1200 
points.  The  DNS  result  shows  a  similar  thin  region  of  reversed  flow  along  the  lower  surface  of  the 
forebody,  near  the  lower  front  comer.  Though  not  shown  in  Figure  19,  the  Cp  distribution  for  the 
DNS  is  similar  to  the  trip  less  DES  prediction  at  Re  =  1  x  105. 

5  Summary  and  Perspectives 

DES  was  applied  to  prediction  of  the  separated  flow  around  an  idealized  jet-fighter  forebody  at  10° 
angle  of  attack.  Influences  of  domain  size,  grid  refinement,  and  turbulence  model  were  investi¬ 
gated  in  cases  in  which  boundary  layer  separation  was  either  laminar,  or  turbulent.  The  initial  and 
boundary  conditions  on  eddy  viscosity  set  the  type  of  boundary-layer  separation,  aimed  at  flows 
above  or  below  the  critical  Reynolds  number  which  controls  a  reversal  of  the  side  force  in  the 
experimental  measurements  of  Polhamus  et  ciL  (1959). 

In  general,  DES  predictions  of  the  super-critical  flow  seem  robust,  tending  towards  experimen¬ 
tal  measurements  with  grid  refinement,  for  example.  The  complex  shedding  process  and  mod¬ 
ulation  in  the  forces  are  consistent  with  circular-cylinder  behavior  and  appear  to  be  represented 
reasonably  adequately,  based  on  the  agreement  between  simulation  and  measurement  of  the  pres¬ 
sure  distribution  and  forces.  In  the  sub-critical  regime,  no  simulation  technique  applied  during  the 
course  of  this  study  (DES/LES,  DNS,  and  preliminary  calculations  using  vortex  methods)  yielded 
a  sign  change  in  the  side  force.  Small  differences  in  the  geometry,  hysteresis,  and  sidewall  effects 
are  three  sources  that  might  explain  the  differences  between  predictions  of  the  sub-critical  flow  and 
experimental  measurements. 
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Aside  from  these  effects,  boundary  layer  treatment  of  the  sub-critical  flow  may  be  simplistic, 
although  fully  “state-of-the-art”.  The  advantage  of  the  tripless  approach  is  that  the  simulation 
parameters  are  unambiguous  (the  results  depend  only  on  the  model,  and  not  on  transition  locations 
chosen  case  by  case).  In  practice,  there  can  be  substantial  regions  of  laminar  flow  and  a  prediction 
of  boundary  layer  transition  is  required.  Details  of  the  processes  of  separation  and  transition  to 
turbulence  as  they  inter-mingle  continue  to  strongly  challenge  current  modeling  approaches. 

In  the  present  application,  the  attached  boundary  layers  lie  entirely  within  the  “RANS  region” 
of  the  DES  solution.  As  the  flow  detaches,  in  the  separating  shear  layers  the  wake  develops  new 
instabilities  that  result  in  the  rapid  growth  of  a  chaotic  and  three-dimensional  wake.  The  lack  of 
eddy  content  in  the  detaching  boundary  layers  represents  a  relatively  small  error  to  the  solutions 
presented  here  (its  length  scale  would  be  several  times  smaller  than  that  of  the  new  instabilities). 
In  other  applications,  e.g.,  flows  with  shallow  separations,  it  will  be  advantageous  and  necessary 
to  seed  the  upstream  flow  with  ‘eddy  content’,  considerably  refine  the  grid  to  support  these  eddies, 
and  initiate  LES  in  the  attached  boundary  layers  prior  to  separation. 

The  current  study  also  provided  an  opportunity  to  assess  the  solver  and  build  confidence  in 
the  application  of  a  second-order  unstructured  method  for  DES  applications.  The  algorithm  was 
sufficiently  accurate  to  capture  the  growth  of  instabilities  in  the  wake  on  both  structured  and  un¬ 
structured  grids  of  reasonable  density.  Streamwise  vortices  were  captured  with  between  five  and 
ten  cells  in  both  cases.  These  features  are  connected  to  aspects  of  the  numerical  method,  such  as 
least-squares  calculations  of  spatial  derivatives  and  the  use  of  non-linear  dissipation  (i.e.,  a  TVD 
type  limiter).  Higher  order  methods  should  be  expected  to  retain  a  given  solution  quality  while 
reducing  the  number  of  cells  required.  However,  obtaining  higher  order  solutions  on  unstructured 
grids  is  significantly  more  challenging  than  on  structured  grids.  Additionally,  higher  order  methods 
would  likely  reduce  the  scalability  of  the  algorithm.  A  side-by-side  comparison  to  a  higher-order 
numerical  method  would  provide  important  estimates  of  potential  benefits  from  a  higher  order  of 
accuracy  (and/or  lower  level  of  numerical  dissipation),  in  turn  allowing  one  to  determine  the  trade¬ 
offs  between  various  approaches.  DES  predictions  on  the  unstructured  grid  showed  the  potential 
benefit  of  the  unstructured  approach  in  placing  points  precisely  where  they  are  needed  -  in  the  near 
wake.  This  study,  however,  is  not  a  comprehensive  testing  of  the  numerics,  since  it  only  included 
experimental  surface  pressures  for  validation.  More  detailed  comparisons,  including  wake  profiles 
and  spectra  are  needed  to  further  assess  the  current  approach. 

The  research  reported  as  part  of  this  undertaking  and  related  efforts  are  increasing  confidence 
in  the  use  of  Computational  Fluid  Dynamics  as  a  more  viable  tool,  capable  of  providing  accurate 
predictions  of  unsteady  flow  in  regimes  which  have  been  either  very  difficult  to  predict  (using 
URANS)  or  pose  too  large  a  computational  cost  (using  LES).  Increasingly,  many  military  aircraft 
are  comprised  of  unusual  configurations,  e.g,  AWACS-like  antennas,  joined  wings,  stealth  shapes, 
etc.  Given  these  disparate  configurations,  there  is  a  pressing  need  for  detailed  information  about 
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their  aerodynamics  that  will  be  continue  to  pose  significant  problems  for  traditional  testing  based 
on  wind-tunnel  and  flight  tests.  CFD  is  positioned  to  provide  detailed  predictions  on  timescales 
that  will  complement  current  efforts  used  for  analysis  and  some  engineering  design. 

Advances  in  DES  currently  depend  both  on  the  physical  modeling,  e.g.,  incorporation  of  rough¬ 
ness  models,  treatment  of  laminar-to-turbulent  transition,  as  well  as  on  numerical  aspects  such  as 
grid  generation,  code  performance  for  unsteady  flow  simulations,  and  issues  related  to  factors  such 
as  numerical  dissipation.  Unstructured  grids  appear  to  hold  great  promise  for  complex  configura¬ 
tions,  a  key  finding  of  the  present  work  was  the  similar,  and  accurate,  prediction  of  the  forces  and 
pressures  acting  on  the  forebody  in  calculations  using  unstructured  grids  as  compared  to  structured- 
grid  results.  This  finding  reinforces  the  importance  of  grid  quality,  e.g.,  low  cell  skewness,  as  the 
important  factor  for  assessing  quality,  rather  than  simply  whether  the  grid  is  unstructured  or  struc¬ 
tured.  For  DES,  guiding  the  grid,  e.g.,  finely  meshing  the  “focus  region”  where  additional  flow 
detail  is  crucial  is  relatively  straightforward  using  unstructured  grids.  Further,  with  the  continued 
development  and  application  of  adaptive  gridding,  the  accuracy  and  fidelity  of  DES  predictions 
should  continue  to  improve. 

References 

[1]  Anderson,  S.B.,  1979,  “Historical  overview  of  stall/spin  characteristics  of  general  aviation 
aircraft”,  J.  Aircraft,  16(7),  pp.  455-461. 

[2]  Batten,  P.,  Goldberg,  U.  and  Chakravarthy,  S.,  2000,  “Sub-grid  turbulence  modeling  for  un¬ 
steady  flow  with  acoustic  resonance”,  AIAA  Paper  00-0473. 

[3]  Bihrle,  W.  Jr.,  Hultberg,  R.S.  &  Mulcay,  W.,  1978,  “Rotary  balance  data  for  a  typical  single 
engine  low-wing  general  aviation  design  for  an  angle-of-attack  range  of  30°  to  90°  ,  NASA 
CR  2972. 

[4]  Constantinescu,  G.S.,  Pacheo,  R.  and  Squires,  K.D.,  2002,  “Detached-Eddy  Simulation  of 
flow  over  a  sphere”,  AIAA  Paper  2002-0425. 

[5]  Durand,  W.F.  (editor),  Aerodynamic  theory:  a  general  review  of  progress,  Springer,  1934-36. 

[6]  Durbin,  P.A.,  1995,  “Separated  flow  computations  using  the  k  -  e  -  v2  model”,  AIAA  J.,  33, 
p.  659. 

[7]  Fritz,  D.,  1995,  “An  investigation  into  the  use  of  two-dimensional  Navier-Stokes  computa¬ 
tions  to  predict  spin  characteristics”,  report  prepared  under  NASA  contract  NAS  1-19672. 

[8]  Fureby,  C.  and  Grinstein,  F.F.,  1999,  “Monotonically  integrated  large  eddy  simulation  of  free 
shear  flows”,  AIAA  J.,  37,  pp.  544-556. 


29 


[9]  Gottlieb,  J  J.  and  Groth,  C.P.T.,  1988,  “Assessment  of  Riemann  Solvers  for  Unsteady  One- 
Dimensional  Inviscid  Flows  of  Perfect  Gases”,  Journal  of  Computational  Physics,  78,  pp. 
437-458. 

[10]  Grismer,  D.S.  &  Jenkins,  J.E.,  1997,  “Critical-state  transients  for  a  rolling  65-degree  delta 
wing”,  J.  Aircraft,  34(3),  pp.  380-386. 

[11]  Hsu,  K.  and  Lee,  S.L.,  1991,  “A  Numerical  Technique  for  Two-Dimensional  Grid  Generation 
with  Grid  Control  at  All  of  the  Boundaries”,  Journal  of  Computational  Physics,  11,  pp.  451- 
469. 

[12]  Hultberg,  R.S.  &  Mulcay,  W.,  1980,  “Rotary  balance  data  for  a  typical  single  engine  general 
aviation  design  for  an  angle  of  attack  range  of  8°  to  90°,  I:  low  wing  model  A”,  NASA  CR 
3100. 

[13]  Jenkins,  J.E.,  Myatt,  J.H.  &  Hanff,  E.S.,  1996,  “Body-axis  rolling  motion  critical  states  of  a 
65-degree  delta  wing”,  J.  Aircraft,  33(2),  pp.  268-278. 

[14]  Jenkins,  J.E.,  1997,  “Nonlinear  aerodynamic  characteristics  of  a  65-degree  delta  wing  in 
rolling  motion:  implications  to  testing  and  flight  mechanics  analysis”,  AIAA  Paper  97-0742. 

[15]  Jobe,  C.E.,  Hsia,  A.H.,  Jenkins,  J.E.  &  Addington,  G.A.,  1996,  “Critical  states  and  flow 
structure  on  a  65-degree  delta  wing”,  J.  Aircraft,  33(2),  pp.  347-352. 

[16]  Karypis,  G.,  Schloegel,  K.  and  Kumar,  V.,  “ParMETIS:  Parallel  Graph  Partitioning  and 
Sparse  Matrix  Ordering  Library  Version  1 .0”,  University  of  Minnesota,  Department  of  Com¬ 
puter  Science,  Minneapolis,  MN  55455,  July  1997. 

[17]  McCormick,  B.W.,  1981,  “Equilibrium  spinning  of  a  typical  single-engine  low-wing  light 
aircraft”,  J.  Aircraft,  18(3). 

[18]  Mittal,  R.  and  Moin,  P.,  1997,  “Suitability  of  Upwind-Biased  Finite  Difference  Schemes  for 
Large-Eddy  Simulation  of  Turbulent  Flows,”  AIAA  J.,  35(8),  pp.  1415-1417. 

[19]  Myatt,  J.H.,  1997,  “Modeling  the  rolling  65-degree  delta  wing  with  critical  state  encounters”, 
AIAA  Paper  97-3646. 

[20]  Nikitin,  N.V.,  Nicoud,  F„  Wasistho,  B.,  Squires,  K.D.  and  Spalart,  P.R.,  2000,  “An  approach 
to  wall  modeling  in  Large-Eddy  Simulation”,  Phys.  Fluids,  12,  pp.  1629-1632. 

[21]  Polhamus,  E.C.,  Geller,  E.W.  and  Grunwald,  K.J.,  1959,  “Pressure  and  Force  Characteris¬ 
tics  of  Noncircular  Cylinders  as  Affected  by  Reynolds  Number  with  a  Method  Included  for 
Determining  the  Potential  Flow  about  Arbitrary  Shapes”,  NASA  TR  R-46. 


30 


[22]  Shur,  M.L.,  Spalart,  P.R.,  Strelets,  M.K.  and  Travin,  A.K.,  1999,  “Detached-Eddy  Simulation 
of  an  Airfoil  at  High  Angle  of  Attack,”  Fourth  International  Symposium  on  Engineering 
Turbulence  Modelling  and  Measurements,  Corsica,  France. 

[23]  Spalart,  P.R.  and  Allmaras,  S.R.,  1994,  “A  One-Equation  Turbulence  Model  for  Aerodynamic 
Flows,”  La  Recherche  Aerospatiale  1,  pp.  5-21. 

[24]  Spalart,  P.R.,  Jou,  W.H.,  Strelets,  M.  and  Allmaras,  S.R.,  1997,  “Comments  on  the  Feasibil¬ 
ity  of  LES  for  Wings,  and  on  a  Hybrid  RANS/LES  Approach,”  First  AFOSR  International 
Conference  on  DNS/LES,  Rouston,  Louisiana,  USA. 

[25]  Spalart,  PR.,  2000,,  “Strategies  for  Turbulence  Modeling  and  Simulations”,  International 
Journal  of  Heat  and  Fluid  Flow,  21,  p.  252-263. 

[26]  Speziale,  C.G.,  1998,  ‘Turbulence  modeling  for  time-dependent  RANS  and  VLES:  a  review”, 
AIAlA  Journal,  36. 

[27]  Steinbrenner,  J.,  Wyman,  N.,  Chawner,  J.,  2000,  “Development  and  Implementation  of  Grid- 
gen’s  Hyperbolic  PDE  and  Extrusion  Methods,”  AIAA  Paper  00-0679. 

[28]  Strang,  W.Z.,  Tomaro,  R.F,  Grismer,  M.J.,  1999,  ‘The  Defining  Methods  of  Cobalt60:  a 
Parallel,  Implicit,  Unstructured  Euler/Navier-Stokes  Flow  Solver,”  AIAA  Paper  99-0786. 

[29]  Strelets,  M.,  2001,  “Detached-Eddy  Simulation  of  Massively  Separated  Flows”,  AIAA  Paper 
2001-0879. 

[30]  Tischler,  M.B.  &  Barlow,  J.B.,  1980,  “Determination  of  the  spin  and  recovery  characteristics 
of  a  general  aviation  design”,  J.  Aircraft,  18,  pp.  238-244. 

[31]  Tobak,  M.  &  Chapman,  G.T.,  1985,  “Nonlinear  problems  in  flight  dynamics  involving  aero¬ 
dynamic  bifurcations”,  AGARD  Symposium  on  Unsteady  Aerodynamics  -  Fundamentals 
and  Applications  to  Aircraft  Dynamics,  AGARD  CP-386,  Paper  25. 

[32]  Travin,  A.,  Shur,  M.,  Strelets,  M.  and  Spalart,  P,  2001,  “Detached-Eddy  Simulations  past  a 
Circular  Cylinder”,  Flow,  Turbulence  and  Combustion,  63,  pp.  293-313. 

[33]  Tromp,  J.C.,  1998,  “Flow  field  simulation  about  a  65-degree  delta  wing  during  constant  roll- 
rate  motions”,  AIAA.  Paper  98-4453. 

[34]  Troung,  K.V.  &  Tobak,  M.,  1990,  “Indicial  response  approach  derived  from  Navier-Stokes 
equations:  Part  1  -  time-invariant  equilibrium  state”,  NASA  TM-102856. 

[35]  Von  Terzi,  D.A.  and  Fasel,  H.F.,  2002,  “A  new  flow  simulation  methodology  applied  to  the 
turbulent  backward-facing  step”,  AIAA  Paper  2 002-0429. 


31 


Figure  1:  (a)  Schematic  representing  coordinate  transformation  from  the  Cartesian  system  to  a 
curvilinear  non-orthogonal  system;  (b)  local  volume  surrounding  the  point  P ,  areas  bounding  the 
control  volume  surrounding  P  given  by  S'7. 
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Figure  2:  Side  view  of  the  rounded-comer  square  (in  a)  and  2:1  ellipse  (in  b ).  Comer  radius  of  the 
rounded-comer  square  is  1/4  of  the  body  width/diameter  D.  Flows  are  from  left  to  right  at  angle 
of  attack  (every  fourth  grid  point  shown  in  the  above  frames). 
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Figure  3:  Side  view  of  the  structured  grid  for  solution  of  the  flow  over  the  rounded-comer  square. 
Comer  radius  is  1/4  of  the  body  width/diameter  D .  Flow  is  from  left  to  right  at  angle  of  attack,  (a) 
entire  computational  domain  for  two-dimensional  simulations;  ( b )  zoomed  view  in  the  vicinity  of 
the  box  showing  clustering  of  grid  points  in  the  wake  region. 
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Figure  4:  Side  view  of  the  structured  grid  for  solution  of  the  flow  over  the  2: 1  ellipse.  Flow  is  from 
left  to  right  at  angle  of  attack,  (a)  initial  domain  for  two-dimensional  simulations;  ( b )  zoomed  view 
in  the  vicinity  of  the  ellipse  showing  clustering  of  grid  points  in  the  wake  region. 
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Figure  5:  Instantaneous  velocity  vectors  colored  by  the  eddy  viscosity  ratio,  Re  —  5  x  104,  10° 
angle-of-attack.  (a)  turbulent  separation  -  inlet  eddy  viscosity  ratio  vtjv  =  3;  ( b )  laminar  sepa¬ 
ration  -  inlet  eddy  viscosity  ratio  vt/v  —  0.  Flow  detachment  near  upper  rear  comer  in  (a),  flow 
detachment  near  upper  front  comer  in  (6). 


Figure  6:  Pressure  coefficient,  Cp,  skin  friction  coefficient,  Cy,  and  turbulence  index  it  around  the 

perimeter  of  the  2:1  ellipse.  Re  =  5  x  104  in  (a),  Re  —  8  x  105  in  (6). - Cv\  e — e  1000  x  C/\ 

- iL.  Angles  are  measured  from  the  fore-symmetry  stagnation  point. 


Figure  7:  Instantaneous  vorticity  contours  in  the  wake  of  the  rounded-comer  square,  Re  —  8  x  10 
turbulent  separation.  Evolution  in  time  through  one  shedding  cycle  evolves  from  (a)  to  (0 


Figure  10:  Three-dimensional  perspective  of  the  structured  grid  around  the  rounded-comer  square. 
Spanwise  dimension  is  six  times  the  forebody  diameter  as  shown  above. 
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Figure  11:  Demarcation  between  RANS  and  LES  regions,  Case  1.  Instantaneous  velocity  vectors 
colored  by  the  eddy  viscosity  ratio  also  shown. 
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Figure  13:  Instantaneous  pathlines  over  the  top  surface  of  the  forebody  colored  by  the  ratio  of  the 
turbulent  eddy  viscosity  to  the  molecular  viscosity  z^/za  The  pathlines  are  shown  from  a  tripless 
solution  in  which  the  upstream  eddy  viscosity  ratio  is  zero  (blue  pathlines  in  the  frame),  non-zero 
values  occur  in  the  separated  region. 
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Figure  14:  Instantaneous  vorticity  contours  in  three  spanwise  planes  along  the  forebody.  Fully- 
turbulent  solutions  at  Re  =  8  x  105,  10°  angle  of  attack,  (a)  Case  1,  (b)  Case  3,  (c)  Case  5,  (d) 
tripless  solution,  simulation  parameters  correspond  to  Case  3. 


Figure  15:  Contours  of  the  instantaneous  vorticity  in  the  forebody  wake  at  Re  —  8  x  105.  (a) 
structured  grid,  Case  3;  ( b )  unstructured  gird,  Case  5. 
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Figure  18:  Temporal  evolution  of  the  streamwise  (Cx)  and  vertical  (Cy)  force  coefficients.  Tripless 
prediction,  simulation  parameters  correspond  to  Case  4. 
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Figure  19:  Pressure  coefficient  distribution  around  the  forebody.  Laminar  separation  cases.  Sym¬ 
bols  are  measurements  from  Polhamus  et  al.  (1959).  - Case  4  (Re  =  1  x  105);  Case  7 

(Re  =  4  x  105). 


